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I. PROTECT OVERVIEW 


I. A. Scientific Objectives In order to understand the physical and chemical processes 

which produce the dust and gas comae in comets and various tenuous planetary and planetary 
satellite (upper) atmospheres through interactions with their particle, field, and radiation 
environs, it is necessary analyze remotely observed and spacecraft data with physically 
meaningful models. With this in mind, we have undertaken a coupled program of theoretical 
modeling and complementary data analysis regarding the distribution of dust and the possible 
production of gas species from dust in comets, and the global distributions of neutral and ionized 
gases in, and escape from, tenuous planetary atmospheres of Io and Europa. The theoretical 
models developed will have further applications to other tenuous atmosphere such as that of 
Pluto or the upper atmospheres of the terrestrial planets. 

I.B. Significance and Relevance This project addresses prime areas of study central to 
current and future NASA missions and for the goals of NASA’s overall program in the study and 
exploration of the solar system and its origin. In the published report of COMPLEX, the 
Committee on Planetary and Lunar Exploration, of the Space Studies Board, National Research 
Council, comets and the Jupiter system were identified in the highest priority group. 

Comets are not only important as members of the solar system and thereby part of 
NASA's mission to explore the solar system, but also in their central role for understanding the 
origin of the solar system and possibly the chemistry of the outer “veneers” of the Earth and 
other terrestrial planets. The study of cometary dust/gas interrelations is important for in situ and 
sample return laboratory measurements in upcoming comet missions. In the near term an 
understanding of sodium dust chemistry obtained through the analysis of remote ground-based 
observations will be relevant in planning for the laboratory studies of harvested cometary dust 
particles in the Stardust Discovery mission, nucleus samples from the planned New Millennium 
Deep Space Four mission, and in situ dust instruments on ROSETTA, an ESA mission to a 
comet with significant NASA participation. 

Io and Europa, their atmospheres, and their interactions with Jupiter’s plasma torus, are 
prominent in current and future NASA science goals and missions. A unique combination of 
modeling tools which have been, and continue to be, developed at the University of Michigan 
have already been applied to understand Galileo particle and fields measurements obtained 
during its encounters with Io and Europa. Both the Galileo mission and the Galileo Europa 
mission and remote measurements with the Hubble Space Telescope will continue to be sources 
of new results, as data already taken is reduced and released, and as upcoming encounters with Io 
and Europa yield new important insight into these atmospheres. NASA is developing plans for 
the Europa Orbiter mission for studying the icy surface and mantle below the icy surface, which 
may hold the clues to the thin extended atmosphere. 

I.C. Strategy and Approach Our overall strategy is and has been to develop and apply 
sophisticated models of the type normally used in purely theoretical studies for the purpose of 
directly comparing models with measurements in order to gain a deeper understanding of the 
physical state and composition of the atmosphere in question. Typical first-principles theoretical 
investigations produce many detailed calculations to explore the type and range of physical 
phenomena that might or should occur, but are often only loosely constrained by observed 
characteristics and are often limited to physically overly simplistic assumptions. On the other 
hand, typical measurements, whether Earth-based and near-Earth based remote observations, or 
in situ spacecraft measurements, are often limited to fairly simple model analyses (scale-height 
distributions, uniform columns, constant temperatures, simple imposed uniform magnetic fields, 
etc.). The Pi’s approach throughout his studies with collaborators dealing both with comets and 
satellite atmospheres has been to use physically meaningful models and compare directly with 
data. For this reason we have developed the core modeling tools, and either analyzed published 
measurements, or established various collaborations with observers. The subjects of study 
covered here are unified by this overarching strategy. 






II. SCIENTIFIC BACKGROUND 

II.A. Scientific Background: Dust and Sodium in Comets The standard dust 
coma treatments originated with the pioneering work of Finson-Probstein (1968a&b) which 
considers the dust to consist of a population of spherical particles having a range of sizes with a 
distribution determined by a power law in particle diameter. The particles were considered to 
have a constant density. p( a) = p i( independent of particle size. The light scattering properties, 
which are needed to describe the observed scattered radiation, and to calculate the radiation 
pressure acceleration on each particle size, were considered to be simply proportional to the 
geometric cross section of the particle implying that the light scattering efficiency, Q s .. h and the 
radiation pressure efficiency, q rp . both equal unity and are independent of particle size. It is clear 
from the work of a number of investigators that these underlying assumptions are highly 
oversimplified for physical dark absorbing dust particles (e.g. see the review by Grim 
&Jessberger 1991). 

The Finson-Probstein approach has been advanced in recent years to include various 
complex physical parameterizations including porous (fractal aggregate) size-dependent dust 
density, realistic radiation and radiation pressure scattering efficiencies for dark dielectrics, 
anisotropic production, and fragmentation (Fulle et al„ 1993 & 1997; Konno et al., 1993; Combi 
1994). The difficulty in interpreting observed monochromatic two-dimensional images of the 
coma and tail is that the parameter space is underconstrained. Even for complicated inversion 
techniques (Fulle et ah), it is necessary to adopt some of the parameterizations, fix some 
parameters, and solve for others. So, one can get a mathematical inversion that is only unique 
within the constraints of the underlying assumptions and parameterizations. 

We (Combi 1994) had found that dust images of the coma provide evidence that the 1-to- 
1 -to- 1 Finson-Probstein type relationship between particle size, velocity, and radiation pressure, 
implicit in most dust coma analyses, does not hold in general. We arrived at a similar and 
independent suggestion as Fulle et ah (1993 & 1997) had: that particle fragmentation within 
and/or just outside the dusty-gas acceleration region (r < 1000km) would produce this effect. 
Analyses of the Giotto Halley Multicolor Camera images independently point to fragmentation 
(Konno et ah 1993; Keller et ah 1990). More recently, Crifo and co-workers (preprint) have 
suggested another possible mechanism which would violate the 1 -to- 1 -to- 1 Finson-Probstein 
relation, that is, if cometary particles come in a range of shapes (i.e., not just spherical, e.g., flat, 
or long and thin), a velocity range for any particular average particle size could also result. It is 
only through continued work on the detailed model analysis of observational data of dust in 
comets, preferably multi-wavelength observations, that further progress can be made. 

Sodium in comets has received much recent attention with the announcement of 
spectacular sodium tail images in comet Hale-Bopp (199502) in the press releases and first 
papers by Cremonese et ah ( 1997) and Wilson et ah (1998). However, sodium tails in comets 
have been reported in observations going back to the early part of the century. Our recent paper 
(Combi et ah 1997) contains a review of many of the older observations. 

II.B. Scientific Background: Tenuous Atmospheres The nature of the atmospheres 
and ionospheres of Jupiter’s natural satellites lo and Europa and their interactions with their 
surrounding radiation, and particles and fields environments is a very active and timely field of 
study. Various kinds of work, depending on different regime-dependent approaches have been 
adopted in recent years, with the hope of understanding the basic global structure of the 
atmospheres, and their interactions w'ith solar radiation and with the Jovian plasma torus 
environment. 

Io’s interaction with Jupiter's corotating plasma torus has been studied for over 25 years. 
Io has a neutral atmosphere which is probably locally thick but rather uneven across its surface. 
(See Lellouch 1996 for an excellent review of pre-1996 literature.) The ultimate source for 
atmospheric gases appears to be the numerous active volcanoes on the surface, moderated by 
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condensation and sublimation from the surface. The energetic particle environment near Io is 
responsible for the balance of the plasma heating, Joule heating, ionization, and surface and 
atmospheric sputtering, and in some form drives the escape of the neutral atmosphere 
(Schneider, Smyth, & McGrath. 1989). 

Our understanding of the interaction of Io with Jupiter’s corotating plasma goes back to 
the discovery of Io-related decametric radio emission discovered by Bigg (1964) and the unipole 
inductor model which explains it (Goldreich & Lynden-Bell 1969). There were a number of 
theoretical studies of this interaction during the immediate post-Voyager era (See Hill et al. 1983 
for an excellent post-Voyager review). Early theoretical work was often done either in the 
context of a “thin” atmosphere [e.g., see Cloutier et al. 1978] indicative of the surface 
temperature (-130 K), or "thick" extended neutral atmosphere (e.g. see Goertz 1980) more 
indicative of volcanic temperatures (-1000 K). Subsequent evidence (again see the review by 
Lellouch 1996) seems to indicate a mixed picture of the global atmosphere, which has a large 
extended corona like a thick atmosphere, but appears to be dominated by local major injection of 
hot (high speed) gas/dust plumes to high altitudes but only near active volcanic vents. Therefore, 
although the atmosphere is probably only locally thick, it still has a large extended neutral corona 
which might provide a sufficient source of impact ionization and photoionization. 

S0 2 appears to be the major primary constituent in lo’s atmosphere (Ballester et al., 1994; 
Lellouch et al., 1992), but observations to date still only provide hemispherically averaged total 
column abundances. The sodium eclipse measurements of Schneider et al. (1991) are the only 
unambiguous information we have about the vertical distribution of neutral gas species in the 
corona of Io from about 1.4 to about 6 R Io . They found a power law with an exponent of a - -3.5 
described this distribution. Smyth & Combi (1997) have shown that this vertical distribution of 
sodium from the Schneider et al. eclipse observations and the more typical emission profiles of 
more spatially extended sodium can together be explained by a single modified incomplete 
cascade sputtering distribution, which results naturally from multiple neutral-neutral and ion- 
neutral collisions between torus ions and atmospheric neutrals. The interpretation of recent 
Hubble observations of various UV emissions of atomic S and O (Ballester et al., 1997) is 
complicated by the fact that the spatial distributions of the emissions are determined by a 
combination of the mass distribution in the atmosphere with excitation by the plasma density and 
temperature environment. Theoretical and empirical models have explored the structure of the 
neutral atmosphere and its behavior to various gas sources and energetic heating and ejection 
processes, such as sublimation, volcanic ejection, surface and atmospheric sputtering, UV, 
plasma, and Joule heating, and LTE and non-LTE IR radiative cooling (Sieveka & Johnson 
1985; McGrath & Johnson 1987; Moreno et al., 1991; Strobel et al., 1994; Pospiezalska & 
Johnson 1995; Wong & Johnson 1995 & 1996). Recent work by Summers & Strobel (1996) 
models the photochemistry and energy balance for lo’s atmosphere accounting for vertical 
transport via the adoption of eddy mixing (and the eddy diffusion coefficient) in a ID-spherical 
model. Wong and Johnson ( 1995 & 1996) have addressed 3D structure as well as horizontal and 
vertical transport explicitly with a hydrodynamic approach. Finally, Marconi et al. (1996) have 
addressed aspects of transition from a lower thick atmosphere to an upper escaping atmosphere 
using a layered approach for a 1 D-spherical, single-species Io-type calculation solving 
hydrodynamic equations at low altitude and performing Direct Simulation Monte Carlo 
calculations at high altitude. Boundary conditions are passed back and forth between the two 
approaches. 

Recent observational studies (Belton et al. 1996; Ballester et al. 1997; Trauger et al. 
1997; Roesler et al. [private communication]) seem to indicate that the atmosphere/ionosphere 
around Io envelopes the satellite at low altitude, but that there are major aspherical anomalies at 
higher altitudes both toward and away from Jupiter, and in the direction of the wake. There may 
be variations with the local magnetic system-III longitude (i.e., the position of the plasma torus). 
Despite all the recent advances in theoretical work, still only a rough sketch of the some aspects 
of the structure of the global atmosphere can be made. It is clear that vertical and horizontal 
transport, day/night and upstream/downstream asymmetries, local peculiarities, and time- 
variations are all important in shaping the global atmosphere. 


4 




There have been a number of 3D numerical studies of the plasma flow past Io (Linker et 
al., 1988, 1989, 1991; Wolf-Gladrow et al., 1987). Because of the complexities in the problem 
and numerical difficulties, they have been limited to restricted portions of the problem, have 
introduced fixes like artificial viscosity, or have been run using less than fully physically realistic 
conditions. Linker et al. (1991) have performed the best example to date of 3D MHD modeling 
of Io. They used the two-step Lax-Wendroff method. These calculations were performed on a 
non-uniform but fixed structured grid, and artificial dissipative terms (artificial viscosity) had to 
be introduced to the basic MHD equations in order to make them computational stable. Our 
contribution in this area is discussed in the next section. 

Although neutral oxygen (Hall et al. 1995) and sodium (Brown & Hill 1996) have only 
recently been detected in a neutral atmosphere around Europa, suggestions that an oxygen 
atmosphere should accumulate above the water-ice surface actually go back 20 years (Yung & 
McElroy 1977; Kumar & Hunten 1982). Further confirmation of the existence of an atmosphere 
was provided recently by detection of an ionosphere in radio occultation experiments from the 
Galileo Orbiter (Kliore et al. 1997). The current working scenario is that water molecules 
breaking up near the surface would, through chemical reactions, produce O, and H,, but the 
lighter H 2 would prefentially escape leaving an 0 2 lower atmosphere, an atomic O upper 
atmosphere/exosphere, and an CP ionosphere. Recent ID-spherical calculations (Ip 1996) suggest 
a possible two component exosphere: a lower thermal population with a small scale height, and 
an upper corona extending a few satellite radii. 

Recent theoretical work on the atmosphere/exosphere/ionosphere of Europa has been 
done using diffusive photochemical and energy-balance models. Being constrained by the 
observations of Hall et al. and Brown & Hill, Ip et al. (1996) have made ID-spherical 
calculations that predict a lower atmosphere with a moderate exponential scale height and an 
extended corona. Sauer et al. (1996) have begun a study of 3D aspects of the interaction of 
Jupiter’s plasma torus with Europa’s atmopshere. They adopt multiple ID hydrostatic-t\ pc 
diffusive profiles for the atmosphere and account for chemistry and energy balance and the 
impacting plasma torus. They solved for the effects of currents and electric fields assuming a 
constant and uniform magnetic field. The recent flybys of Europa by the Galileo orbiter indicate 
the presence of an ionosphere (Kliore et al. 1997), roughly consistent with expectations from the 
O detection, and a disturbed local magnetic field (Kivelson et al. 1997) which could bo 
accounted for by an intrinsic field or possibly an induced field (Kivelson, private 
communication). Our new studies of the effects of ionospheric mass-loading on the flow and 
local magnetic field are discussed in the next section. 

III. PROGRESS AND CURRENT WORK STATUS 

III.A. Progress of Work: Dust Coma Studies During the last few months of the 
current grant, our plan was to complete the analysis of the comet Austin (1989 Cl) dust coma 
images in collaboration with Dr. David Schleicher. For this work we are performing a similar 
dust coma analysis as was done for the comet P/Halley images as described in the paper by 
Combi (1994). Observations of the coma of comet Austin indicated a fairly steady and uniform 
production of gas and dust (Schleicher, private communication), not complicated by large 
amplitude time-dependent variations or aspherical production. It provides a good testbed for 
studying the role for dust velocity dispersion (for each particle size) that could result from 
particle fragmentation. We have contour plots of one of the images, however, Dr. David Osip 
who was a postdoc with Dr. Schleicher during the observation and reduction of comet Austin 
data, was to provide us with the images in digital form. In anticipation of receiving the digital 
data and using the contour plots, the modeling background work was begun and preliminary 
model calculations were performed and compared with the observed coma. 

For Halley, Combi ( 1994) found that 3 images of the dust coma could not be adequately 
reproduced with a straightforward Finson-Probstein type formulation. A model which produced 
dust particles that fragmented into smaller particles after they reached terminal velocity in the 
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dust-gas drag region could reproduce the images. In keeping with the Monte Carlo nature of the 
calculation a “heuristic” dust fragmentation parameter was constructed in the form 

df = d 0 R a 

where d G is the radius of a dust particle before fragmentation, R is a random number on the 
interval of 0 to 1, and a is an adjustable postive parameter. There is no fragmentation for a value 
of a = 0, and much fragmentation for a value as large as 4. This formulation has the interesting 
property of conserving the form of the size distribution which varies from d' V2 . to d' ? 4 . Combi 
(1994) had found that with a value of a = 2 three images of the dust coma taken at heliocentric 
distances of 1.82, 1.45 and 1.08 AU before perihelion could all be well reproduced by the model. 
A value of p = 2 means that each original particle fragments on the average into 27 particles. 
This is consistent with the analysis of the Giotto camera images of the dust coma very near to 
Halley’s nucleus by Thomas and Keller ( 1990), who found that each original dust particle had to 
fragment into at least 2.7 particles as a lower limit. 

For comet Austin the first sets of model calculations indicated that the power law 
exponent of the particle size distribution had to be changed. This is actually consistent with the 
results of Fulle et al. (1993) for comet Austin. It is clear from the preliminary results that a 
model with a moderate amount of fragmentation will be required (a ~ 2) but without a careful 
comparison with the actual digital image data it is impossible to reach an clear conclusions at this 
point. This work will be continued and completed in the coming months after the current grant 
period, but the preliminary comparison is illustrated in Figures 1 and 2. 

In addition to the apparent direct production of gaseous sodium (Na) from the nucleus 
(Combi et al. 1996, 1997; Cremonese et al. 1997; Brown et al. 1998), there have been reports of 
observations of comet Hale-Bopp giving convincing evidence for a possible production of 
gaseous Na from dust, both in the coma of the comet (Brown et al. 1998; Rauer et al. 1998; 
Arpigny et al. 1998), and far down the dust tail (Fitzsimmons et al. 1997; Wilson et al. 1998). 
Figures 3 a&b show two images of the cometary Na tails taken by Wilson et al. in March 1997 
and by Cremonese et al. in April. Both the brightness down the tail and the tail position angle 
with respect to the sunward direction from the April image has been shown with some simple 
model calculations by Cremonese et al. to be consistent with a nucleus or near-nucleus source, 
similar to the direct nucleus source component seen in comets observations of comets Bennett. 
Kohoutek and Halley (Combi et al. 1997). However, Wilson et al. conclude from their March 
image that neither the brightness distribution nor the tail position angle are consistent with a 
direct nucleus or near nucleus source. Some sodium must be produced by an extended source 
out in the dust tail. Based on long-slit spectroscopy of the coma region (distances < 2 x 10 ? km) 
Brown et al. (1998) report roughly equal contributions to the total gaseous Na production rate in 
comet Hale-Bopp coming from the direct nucleus source and the extended dust source. This 
raises interesting questions regarding the processing of dust particles once liberated from the 
nucleus, both over short periods of time in the coma and long periods of time required to travel 
large distances down the tail and to fill interplanetary space. Is observed cometary Na eroded 
only from the surfaces of typical dust particles? Or, does the gaseous Na originate from a 
population of dust particles, like the mixed CHON-refractory particles that might completely 
disintegrate? The careful detailed modeling of remote ground-based observations of gaseous Na 
in the comas and tails of comets will provide important results in this area. 

Spatial profiles of Na in three comets (Bennett 1969 Yl, Kohoutek 1973 F, and 
lP/Halley) were determined by Combi et al. (1997) from long-slit spectra (2D spatial/spectral 
images) with the slit aligned both parallel and perpendicular to the comet-sun line, and displaced 
antisunward. Like the various reported Hale-Bopp measurements they found evidence for a 
source of gaseous Na which originated at the nucleus (or a very short-lived parent species). 
Unlike Hale-Bopp which seems to show an extended source spatially coexistent with the dust in 
the sunward coma and tail, we found a large and highly temporally variable extended source on 
the tailward side of the nucleus, that bore a remarkably similar spatial morphology to that of 
H,0 + , but no similarity to the dust continuum profiles determined from the same data. This 
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Figure 1. Fragmentation Models for Comet Austin. Shown from the upper left are four models for the 
comet Austin observation geometry in mid-May 1990 with four values of the fragmentation parameter, a. As in the 
case of Halley, clearly the simple model without fragmentation (a=0) is not correct, hut at this point it not clear 
which model is best when comparing with the contour plot of the observation in Figure 2. 

tailward extended source varied by factors of a few from day-to-day in two of these comets, and 
was most noticeable in observations of these “old" comets at heliocentric distances < 0.9 AU. 
All three of these comets had similar gas production rates that were more than an order of 
magnitude smaller than Hale-Bopp at comparable heliocentric distances. Clearly an important 
aspect of understanding cometarv sodium will be differences from comet to comet. 

Combi et al. showed that a direct nucleus source model could describe the sunward 
profiles and the inner part of the antisunward profiles in the three “old” comets. It accounted for 
direct production of Na from the nucleus, the heliocentric velocity dependences of the radiation 
pressure acceleration and the g-factor (from the solar Fraunhofer lines), collisional entrainment 
within the inner part of the expanding water coma, and photoionization. We showed that (1 ) 
collisional entrainment could explain the sunward extension of the profile, (2) an extended 
source was required to explain most of the tailward spatial profiles in all three comets, (3) the 
photoionization lifetime as computed from theoretical cross sections, which is 3 times longer 
than that computed from old experimental cross sections (Huebner, Keady and Lyon 1992). 
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seems to be correct, and (4) that the nucleus (or near-nucleus) source of gaseous Na seen in the 
coma was about -1% of that expected from solar abundances. Cremonese et al. (1997) 
independently reached a similar conclusion about the Na lifetime and the production rate in 
comet Hale-Bopp. All the observational and model analysis results of Combi et al. (1997) had 
already been presented at the Asteroids. Comets, Meteors conference by Combi et al. 1996, 9 
months before any of the observations of Na in comet Hale-Bopp. 



Figure 2. Contour Plot of Dust Coma Image of Comet Austin. Shown is a contour plot of the dust 
coma of comet Austin taken with a 4845 A filter. The square field of view is roughly comparable to the models 
shown in Figure 1 . 


Brown et al. (1998) have compared a Monte Carlo model to their measured spatial and 
line-of-sight Doppler velocity profiles of Na in comet Hale-Bopp along a long slit aligned with 
the comet-sun line. Their model includes a calculation for an extended dust source and a nucleus 
source. Their analysis notes and indirectly accounts for the collisional coupling between Na 
atoms and molecules in the inner coma by avoiding that portion of the calculation. None of the 
published models so far (Combi et al. 1997; Cremonese et al. 1997; Brown et al. 1998; Kawakita 
and Fujii 1998) can address all combinations of effects which are actually important to sort out 
all the important physical effects at the same time, in all the data, including the orbital dynamics 
effects and extended production from the dust coma and tail. 

Work has begun to incorporate an extended source of Na from dust in our full 3D dust 
model and to make a first study. The shell of the dust coma model is actually the same code as 
was developed years ago for our hydrogen coma model (Combi & Smyth 1988a&b). A 
calculation of an extended Na from dust model then is just a merger of the core source region 
parts of the H coma model and the dust coma model. In this calculation the time-history of a Na 
atom follows the same sequence as now done for H atoms. The calculation accounts for the 
extended production of Na atoms from a dust coma/tail distribution, the time-dependent variation 
of the dust/gas production rates, the heliocentric distance dependence of the Na lifetime, the 
heliocentric distance and velocity dependence of the Na g-factor and solar radiation pressure 
acceleration. These results can then be added as necessary to a direct nucleus component for 
comparison with data which of course includes Na from both sources. 
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March 17, 1998 


April 17, 1998 



Figure 3. Images of the Sodium Tail in Comet Hale-Bopp. At the left is a contour plot of the sodium 
image from Wilson et al. (1998) alter complete subtraction of the underlying dust superimposed over the false color 
versions of the H 2 CT and dust-continuum tails. At the right is the narrow-band interference filter image obtained by 
Cremonese et al. (1997) which includes a combination of sodium (the long thin feature running from the head 
toward the uppermost right) plus the underlying dust tail (the diffuse structure running toward the right). Clearly the 
long, thin and bright sodium tail seen by Cremonese et al. in April was not apparent a month earlier. 

Figure 4 shows our prelimary model calculation for the Na distribution for comet Hale- 
Bopp as produced by the dust in the coma and tail, corresponding to the March 17, 1997, data of 
Wilson et al. (1998). Figure 5 shows a direct nucleus source model for the April (Cremonese) 
observations. The direct nucleus model does produce a long straight Na tail as described by 
Cremonese et al. (1997) for April. (Note the spatial scale of the model is much larger than the 
image. The observations only cover the inner straight portion of the Na tail which is shown in the 
model to have some curvature at larger distances.) The model reproduces both the variation of 
brightness down the tail and the position angle of the tail in a single self-contained and self- 
consistent calculation. Cremonese et al. had presented some calculations for the average velocity 
and brightness down the tail using a simplified calculation and compared it with their spectral 
observations. Separately, they presented some syndynamic calculations of the apparent position 
angle of the sodium tail in their image data, given some representative values of the p parameter 
(the ratio of the radiation pressure acceleration to the sun’s graviational acceleration). Our 
model combines all effects realistically in a single Monte Carlo numerical simulation As it 
should be, our nucleus-centered model is consistent with the simple model calculation of 
Cremonese et al. (1997) in terms of the brightness and velocity variations down the tail, as well 
with as the position angle of the Min. 
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Figure 4. Model Calculation for a Dust Source of Gaseous Sodium. Shown are contours ol 
brightness (in kilorayleighs) for a model calculation where sodium atoms are produced from dust with a production 
rate proportional to the surface area of the dust particles. The sun is to the left. 

The direct nucleus source model for the March (Wilson) observation, not shown here, 
looks very similar to that for April (although overall not as bright as April for the same 
production rate). Clearly it does not look anything like the March observation. On the other 
hand, the basic behavior in the March image is reasonably well reproduced by the model 
produced with the dust source. No effort has yet been made to carefully reproduce the brightness 
distribution or simultaneously reproduce the dust tail image by adjusting the dust model 
parameters. This work is left for after the current grant period. The dust model includes a 
nominal dust size power law distribution (28 particle sizes distributed logarithmically from 0. 1 
pm to 1 cm). Dust terminal velocities are calculated using Sekanina’s approximate formula (see 
Sekanina 1981 and Sekanina and Larson 1984). The model reproduces the approximate position 
angle of the Na tail (about mid-way between the ion and dust tails) and the small peak in 
brightness (750 Rayleighs) located at about 2 million km down the tail as seen in the 
observations of Wilson et al. (1998). This peak results from sodium atoms produced in the 
vicinity of the nucleus and pushed to large enough antisunward velocities so that they appear tin 
velocity space) to go first through the center of the solar Fraunhofer line absorpton and over to 
the red side wing. This produces a large increase in the effective g-factor between the nucleus 
and this location. (See also the later discussion in connection with Figure 7.) 
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Figure 5. Model Calculation for a Direct Nucleus Source of Gaseous Sodium for April 24. 
1997. Shown is a model calculation for a direct nucleus source for sodium under the influence of collisional 
entrainment in the outflowing water coma, the heliocentric velocity dependent radiation pressure acceleration and g- 
factor, and 3D time-dependent solar orbital dynamics. The linear dimension is as in Figure 4: 100 units = 2 million 
km. The sun is to the left. The map is projected for the Earth-comet-Sun geometry of the date. 

Figure 6 shows the result of a calculation of the dust source model corresponding to the 
April postperihelion image of Cremonese. The sodium produced from dust in the inner coma 
does produce a long straight tail similar to that observed, but the brightness profile close to the 
nucleus for this extended source does not reproduce the observations. However, based on the 
Cremonese et al. image alone, it is not yet possible to put a limit on the contribution to the total 
sodium product in the dust source compared with the direct nucleus source because the dust 
contribution was not subtracted. 

Recently, Kawakita and Fujii (1998) have presented model calculations for the March 
(Wilson) and April (Cremonese) images using the same Monte Carlo procedures as had been 
applied by us (Combi et al. 1997) for the sodium distribution in the inner coma of the three “old" 
comets. This calculation includes the heliocentric velocity dependent radiation pressure 
acceleration and g-factor for sodium atoms and collisional entrainment of sodium atoms in the 
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Dust Source for Na in Hole— Bopp: April 24, 1997 (kilo ray lefqhs) 



Figure 6. Model Calculation for a Dust Source of Gaseous Sodium for April 24, 1997. Shown 
is a model calculation for a dust source of gaseous sodium corresponding to the image of Cremonese et al. ( 1997) 
who report a long straight tail at this time. 

outflowing water-dominated coma, but does not include the full 3D sun-centered orbital 
dynamics calculations of our new direct nucleus model, nor does it address the dust coma/tail 
source for sodium. Kawakita and Fujii conclude (based on these rather incomplete modeling 
studies) that the major difference between the March and April images can be explained simply 
by the difference in the radiation pressure acceleration environment between the two periods. 
There is some validity to their suspicion in this regard, but it does not explain the whole story. 
Figure 7 shows a plot from Combi et al. (1997) showing the variation of the radiation pressure 
acceleration on sodium atoms with heliocentric velocity. Note that the g-factor, not surprisingly, 
also has a similar shape. The heliocentric velocity of comet Hale-Bopp for the March and April 
time periods are also noted in the figure. Note the arrows anchored at the positions of the 
comet’s rest frame for each date show the direction (red shift) of liberated sodium atoms in the 
coma and tail. It is true that one would expect the tail to be longer and the brightness to be 
higher in April than March as suggested by Kawakita and Fujii. However, their calculation does 
not (and cannot) address the difference in position angle with respect to the sunward direction 
between the two time periods. The difference in position angle is not explained simply by the 
difference in heliocentric velocity. So, although they may be correct in stating that the long 
straight Na tail produced by the direct nucleus source should be weaker in March than in April, 
the appearance of the overall tail in March cannot be explained simply by this effect alone. It is 
clear based on our preliminary calculations, which are consistent with the suggestions of Wilson 
et al. (1998). that a direct nucleus source alone cannot explain the position angle of the tail. Our 
preliminary dust coma/tail source model (Figure 4) clearly shows that a dust coma/tail source 
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probably explains most of sodium that is seen in the extended tail images of Wilson et al. in 
March. 


Radiation Pressure 
Acceleration 



Figure 7. Heliocentric Velocity Dependence of the Radiation Presssure Acceleration and 
Sodium G-factor Shown is a figure from Combi et al. (1997) of the variation of the radiation pressure 
acceleration (above) and fluorescence rate (below) with an atom's heliocentric velocity. The heliocentric velocity of 
an atom in comet Hale-Bopp's rest frame for mid-March (Wilson et al.) and mid-April (Cremonese et al. image) 
1997 are shown. 
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Finally, the PI was invited to give a talk on Dust-Gas Interrelations at the First 
International Conference on Comet Hale-Bopp which was held on February 1998. The topic 
involved dust coma subjects appropriate to this project and gas coma subjects appropriate for 
another. Work in the months preceding the talk involved reviewing dust and gas issues 
pertaining to comets in general and Hale-Bopp in particular, such as extended source of gas from 
dust, dust fading, and dust fragmentation. A paper “Dust-Gas Interrelations in Comets: 
Observations and Theory” (Combi et al. 1998b) was submitted for the special refereed 
proceedings volume of “Earth, Moon, and Planets As of this writing the paper has been revised 
and is about to be resubmitted. From the language in the review and interchanges with the editor 
we expect this version will be accepted in the next few weeks. 

III.B. Progress of Work: Tenuous Atmosphere Studies Two models and 
modeling efforts serve as the bases of the work described in the tenuous atmosphere portion of 
this proposal. One is a major modeling effort directed by our colleague. Dr. Tamas Gombosi, 
involving several complementary subsets and collaborative efforts with other independent 
scientists at Michigan (like the PI and Co-I of this proposal), to develop robust and efficient 3D 
MHD and hydrodynamic computer codes for application to various problems in Earth and Solar 
System space science. These include, for example, the Earth’s magnetosphere, the heliosphere 
(Linde et al. 1998), the comet solar-wind interaction with ion-chemistry (Gombosi et al. 1996; 
Haberli et al. 1997), the dusty-gas comet coma (Combi et al. 1998b), the comet-like interaction 
of the solar wind with Mars and Venus (Bauske et al. 1998), the magnetospheres of the outer 
planets, and the interactions of planetary satellite (lo, Europa, and Titan) atmospheres with the 
local corotating planetary plasma (Combi et al. 1998a; Kabin et al. 1997abc, & 1998a). 

3D Multiscale MHD Model with Mass-Loading The first results for applying the 3D multiscale 
ideal MHD model for the mass-loaded How of Jupiter’s corotating magnetospheric plasma past 
Io are described in our recently published paper (Combi et al. 1998a). A copy of this paper is 
included as part of this report. Preliminary versions of these results were presented last year by 
Combi et al. (1997) and Kabin et al. (1997a). The model is able to consider simultaneously 
physically realistic conditions for ion mass loading, ion-neutral drag, and intrinsic magnetic field 
in a full global calculation without imposing artificial dissipation. Io’s extended neutral 
atmosphere loads the corotating plasma torus flow with mass, momentum, and energy. The 
MHD equations for mass, momentum, energy and B-field are solved on an adaptively refined 
unstructured Cartesian grid using a multiscale adaptive upwind scheme for MHD (MAUS-MHD) 
with a novel scheme for dealing in a numerically robust way with the difficulties normally 
encountered in implementing the div B = 0 condition (Powell 1994). For the work described our 
paper we explored a range of models without an intrinsic magnetic field for Io. Figure 8 shows a 
comparison of our results with particle and field measurements made during the December 7, 
1995, flyby of Io, as published by the Galileo Orbiter experiment teams. For two extreme cases 
of lower boundary conditions at Io, our model can quantitatively explain the variation of density 
along the spacecraft trajectory and can reproduce the general appearance of the variations of 
magnetic field and ion pressure and temperature. The two cases correspond to boundary 
conditions at an ionosphere boundary of 150 km above the surface of Io for the plasma density 
and the B-field which are fixed (an absorbing boundary, allowing plasma flow and magnetic flux 
to pass through) and reflective (a robust non-penetrating and perfectly conducting ionosphere). 
The effect of a realistic ionosphere is somewhere between these extremes. The net fresh ion 
mass-loading rates are in the range of -300-650 kg s' 1 , and equivalent charge exchange mass- 
loading rates are in the range -540- 1 1 50 kg s' 1 in the vicinity of Io. 

First calculations have also been performed for the interaction of Europa’s atmosphere 
with Jupiter’s plasma torus ( Kabin et al. 1997b). A paper has been submitted to JGR (Kabin et al. 
1998b) and is currently being revised. A copy of the submitted version is included as part of this 
report. As a starting point we adopted the simple atmosphere derived from the Galileo radio 
occulation results of Europa’s ionosphere by Kliore et al. (1997). However, it is necessary to 
point out that there are large uncertainties about the mean of the suggested ionosphere, and the 
derived neutral atmosphere from Kliore et al. A given neutral atmospheric surface density and 
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scale height, combined with a net ionization rate provides ion mass, momentum and energy 
loading rates, by analogy with our Io results, and serve as a source term for the MHD equations. 
We used the vacuum superposition approach of Kivelson et al. (1997) as a starting point for 
estimating an intrinsic magnetic field for Europa and then iteratively constrained our model 
parameters by comparison with the Galileo magnetometer measurements (Kivelson et al. 1997) 
which we have obtained in numerical form from Dr. Kivelson and co-workers at UCLA. Not 
surprisingly the addition of mass loading, decreases the dipole moment required for an intrinsic 
field. So far our best results have been obtained with a model which has an ion mass loading rate 
to 6 x 10 24 s '. This is below the rate suggested from the nominal mean atmosphere of Kliore et 
al., however, the difference is not unreasonable given the uncertainties in the scale height from 6 
occulations, the values of the appropriate ionization and recombinations frequencies (factors of a 
few), and the simplicity of the spatial nature of our source rates. For our best model the spatial 
extent and magnitude of the enhancement in plasma density near Europa, compared with the 
nearby ambient plasma torus density, is very similar to those shown by Gurnett et al. (1998) for 
electron density as derived from the Galileo plasma wave measurements. The preliminary 
results are shown in Figure 9. 



Figure 8. Comparison of MHD Mass-Loading Models with Galileo Particle and Field 
Measurements. The magnetic field measurements from Kivelson et al. (1996) are in (a), and the plasma density, 
pressure and temperature measurements from Frank et al. (1996) are in (h), (c), and (d), respectively. The solid 
model lines are for reflective boundary conditions and the dashed model lines are for fixed boundary conditions at 
an altitude of 150 km from Io by Combi et al. (1998a). 

Kinetic Ion-Neutral DSMC in 3D ( KIND_3D) The new model undertaking of this project is a 3D 
neutral and ion multi-species kinetic model, embedded in the multiscale Cartesian grid system 
structure of the 3D MHD code. The KIND_3D is based on DSMC techniques taken from an 
extensive repertoire of dilute gas dynamics methods in the computational fluid dynamics (CFD) 
literature, and already applied to cometary and Earth polar wind applications by our group. For 
the original CFD applications see, for example, the monograph by Bird (1994). The method 
enables a kinetic calculation of neutral and ion gas species to be made, either separately, 
iteratively, or ultimately together, in a satellite or tenuous planetary atmosphere, ionosphere, and 
exosphere/corona. Results from either MHD and DSMC models can then easily serve as inputs 
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to the other, or ultimately the models can be run intimately. Work in Michigan on kinetic 
modeling of neutral and ionized gases using particle simulation codes began in the last few years 
in applications for multispecies 2D axisymmetric modeling of the non-equilibrium flow of the 
neutral coma (Combi 1996), with non-LTE radiational cooling and UV heating, and for 
multispecies ID-spherical modeling of the Earth’s polar wind flow (Miller & Combi 1994; 
Miller et al. 1995). A description ofthe KIND_3D model for Io was presented at the 1998 Spring 
AGU meeting in Boston (Bauske & Combi 1998a). First calculations were presented at the 1998 
Fall AGU meeting in San Francisco (Bauske & Combi 1998b). For these results we ran a single 
species neutral model for an atmosphere produced by a combination of sublimation and two 
volcanoes. A major modeling paper, documenting KIND_3D and presenting first applications to 
Io is being prepared for submission during the early part of the follow-on grant which has 
already been approved for funding by the Planetary Atmospheres. 

The Direct Simulation Monte Carlo (DSMC) method has a long history (Bird 1994) of 
application in areas as diverse as study of hydrodynamic shock structures, calculation of low- 
density, hypersonic flows around spacecraft, and simulations of plasma reactors for 
microelectronics manufacturing. Originally, the method was used for simulations in the (neutral 
gas) transition regime, where the mean free paths of particles are too large for traditional 
continuum CFD methods to be applicable. Because collisions are important, free molecular 
simulations are also not appropriate. Therefore, individual particles are simulated as they move 
around within a grid, colliding with other particles and with solid objects. Macroscopic 
properties, such as density and temperature are computed by appropriate averaging of particle 
masses, positions, velocities, and internal energies. Surface properties are calculated from the 
momentum and energy exchanges with the surface, also allowing for chemical reactions. In the 
collisionally thick limit an appropriate DSMC calculation asymptotically approaches the same 
results as a hydrodynamic calculations such as the Euler or Navier-Stokes formulations. 


Reduced mass-loading and d pole Reduced mass4oadlng and dipole 
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Figure 9. MHD Model with Mass Loading for Galileo E4 Flyby of Europa. At the left are shown 
the 3 components of magnetic field compared with the Galileo magnetometer measurements. At the right are shown 
plots of plasma density, velocity, and temperature. The model has fixed boundary conditions at the surface, mass- 
loading, ion-neutral charge exchange, recombination, and an intrinsic magnetic dipole field with components (0, 7 1 , 
-64) R e 3 nT. Results are from Kahin et al. ( 1997b). 

Many applications today are in dilute, but collisionally dominated particle regimes where 
the particles have nearly thermal collision energies and the ionization ratio is very low. In 
particle oriented plasma physics the Particle In Cell (PIC) method is used, usually in order to 
study collision free plasmas (for a review see Winske & Omidi, 1996) and therefore particle 
collisions were included only in special cases (see, e.g., Wang et al., 1996 and references 
therein). With the hybrid model we explore a new way to bridge the gap between the collision 
dominated and the collision free regime in partially ionized plasmas and we also extend the 
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range of the DSMC method towards the higher particle energies usually encountered in planetary 
ionospheric and magnetospheric plasma environments (hundreds of eV/arau). 

We divide the spatial region simulated by the model into ‘inner’ and ‘outer’ domains. 
Inner domains may contain a whole planetary obstacle and a fraction of its surrounding 
atmosphere / exosphere / plasma environment. Both types of domains are simulated using the 
MHD approach in order to obtain the magnetic field. In inner domains we embed the 
multispecies KIND_3D model which allows detailed simulations of the chemistry and of surface 
interactions. Each domain is adaptively refined and therefore contains multiple cells in various 
sizes to match requirements specific to the domain. The MHD parts, e.g., may have high 
resolution at a planetary bow shock, while the DSMC parts are able to simulate the outflow of 
hot gases due to volcanic activity from specific surface regions above Io. At domain boundaries 
the MHD flow entering a DSMC domain is converted to entering particles by sampling from the 
particle distribution. 

Lutisan (1995) compared several variants of the DSMC method and found the No Time 
Counter technique (NTC) of Bird (1994) and the Null-Collision technique (Koura, 1986, 1989. 
1998) to be among the fastest and most reliable ones. We chose the NTC method for its 
simplicity and its closeness to classic kinetic theory. In contrast to Bird, we allow a particle to 
undergo several changes of its internal energy per time step, thus allowing for fewer particles per 
cell and larger time steps where necessary. 

Apart from higher energies than in usual DSMC calculations, we also face the need to 
simulate regions with orders of magnitude differences in density. The model chemistry will 
contain a variety of different species including trace species. Our basic ideas to approach these 
problems numerically are the use of weighted species groups and of a regular spaced subgrid 
within each cell of the adaptive grid structure. Usually every particle in a simulation cell 
represents a fixed number of real particles. This leads to problems if one species is 
underrepresented since a large number of samples is necessary to get reasonable statistics for this 
species. A weighting method artificially changes the species value in order to allow for more 
simulation particles. In collisions of species with different weights additional means must be 
undertaken to conserve the overall momentum and energy balance at least on average. For binary 
Coulomb collisions and a completely different particle selection method Takizuka & Abe ( 1977) 
suggested an algorithm which was later improved by Miller & Combi (1994). 

KIND_3D has added functionality beyond the basic algorithms suggested by Bird ( 1994) 
to account for chemical reactions and collisions with large inelastic cross sections. We developed 
a corresponding weighting scheme which allows for trace species. Basically the number of 
selections for collisions is now recalculated every time a species type changing collision 
(reaction, ionization, etc.) occurs and the collision iteration proceeds as long as there are 
selections to work on. Our NTC variant employs a data structure which allows fast addition and 
deletion of particles without need to loop over all particles for reindexing after such collisions. 
The weighting scheme determines the number of selections and the acceptance of the collisional 
momentum, energy, and type change for the higher weighted particle. Weighting is also used to 
keep the numbers of simulation particles in different density regimes approximately equal and at 
a computationally feasible number. Care is taken of particles moving from one cell to another 
with a different cell weight. Particles with higher weights are replicated, and particles with lower 
weights are deleted with a probability set by the ratio of the particle weight to the cell weight. 

Now we describe the current status of the model, as of this writing. A single processor 
version of the hybrid model is specified for Io. Particle species include O, S, 0 2 , S 2 , SO, SO ; , 
SO,, and the corresponding ions. We account for 0 + ( 4 S), 0 + ( 2 D), 0*( 2 P), S + ( 4 S), S + ( 2 D), S + ( : P) as 
separate species for photoionization and take care of the excess energies which contribute to the 
electron energies for ionization, to the particle energies for dissociation, and to UV heating. 
Energy dependent solar photo rates and the solar photon fluxes are taken from Huebner et al. 
(1992). Particle trajectories are calculated either assuming the adiabatic particle motion in the 
corotating magnetic field of Jupiter— the gyrotropic approximation— (Northrop & Hill, 1983; 
Northrop & Birmingham, 1982) or by direct numerical integration of the momentum equations. 
The plasma is assumed to be quasi neutral and as a first approximation the electron temperature 
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given by the scaled MHD plasma temperature. At a later stage we will calculate electrons 
independently, for example, using the method of Harvey & Gallis (1996). 

Insufficient knowledge about the energy dependence of collision cross sections is a major 
problem in space physics. KIND_3D includes the variable soft and hard sphere molecule models 
(Bird 1981; Koura & Matsumoto 1991, 1992) for elastic cross sections. For small particle 
energies Bird (1994) suggests schemes to convert temperature dependent reaction rate 
coefficients to steric factors. We are extending these and the elastic cross sections towards higher 
energies with theoretical approximations from Johnson (1990) and will make use of accurate 
calculations from the astrophysical community (e.g., the Opacity project, Cunto et al. 1993) to 
get a fast and reliable chemical reaction network for high collision energies. An up-to-date 
atomic/molecular database will be maintained. 

In DSMC methods dissociation and vibrational excitation during molecular collisions are 
coupled. Whenever a redistribution of vibrational energy during a collision occurs, there is some 
probability for dissociation if the available energy exceeds the dissociation threshold. 
Calculations of these effects are usually based on the statistical collision model of Borgnakke & 
Larsen (1975) and its improvements made by a number of successors. As a first step, we use the 
method suggested by Bird (1994) and Carlson & Bird (1995) for molecular collisions. Here 
discrete vibrational energy levels are given by a harmonic oscillator potential. The energy levels 
are extrapolated beyond the dissociation level and a probability is calculated in order to decide 
wether the molecule dissociates or not. A new scheme was suggested by Lord (1998) which does 
not need the extrapolation and uses a more appropriate anharmonic oscillator potential. We will 
implement it in a second step. Cooling, or radiation escape losses are the other half of accounting 
for internal rotational energy transfer in collisional processes. Strobel et al. (1994) accounted for 
non-LTE radiative escape in the context of their ID diffusive photochemical model. Non-LTE 
radiative escape is handled in DSMC as a by-product of the accounting of internal energy, as 
described for water IR cooling in a cometary coma DSMC model (Combi 1996). 

With the DSMC method a lower boundary can be adopted either from somewhat below a 
defined exobase at some altitude above the surface, or at the surface (as is the case possibly at 
locations on the night side of Io and/or far away from volcanoes, or for Europa). Particle surface 
interactions are currently restricted to elastic collisions or absorption. One of the next steps 
among others will be the implementation of the Cercignani-Lampis surface interaction model 
(Lord, 1991). Depending on the obstacle chosen for the MHD part, the lower boundary is not 
always equal to the planetary surface. In this case incoming particles will be sampled from 
exospheric fluxes which are specified depending on surface conditions. 

KIND_3D provides a self-consistent approach to modeling chemistry, radiative effects, 
and kinetic effects important in the interaction of the extended neutral atmosphere and 
ionosphere with the outer impinging magnetized plasma flow. The advantage of having KJND- 
3D embedded in the MHD multiscale adaptive grid means that both parts of the problem can be 
solved using the most physically reasonable but yet still computationally efficient methods. 
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Abstract. The three-dimensional interaction of Europa with the Jovian 
magnetosphere is modeled using a complete set of ideal magnetohydrodynamic (MHD) 
equations. The model accounts for exospheric mass-loading, ion-neutral charge exchange, 
recombination and a possible intrinsic dipole magnetic field of Europa. The single fluid 
MHD equations are solved using a modern finite-volume higher-order Godunov-tvpe 
method on an adaptively refined unstructured grid, which allows a detailed modeling of 
the region near Europa, while still resolving both the upstream region and the satellite's 
wake. The magnetic field and plasma density measured during Galileo’s E4 flyby of 19 
December 1996 are reproduced reasonably well in the simulation. The general structure 
of the plasma flow past Europa and bounds on the total mass-loading and ion-neutral 
charge exchange rates are discussed. 
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Introduction 

Jupiter’s satellite Europa is located deep inside Jovian magnetosphere on the 
outskirts of the Io plasma torus. Europa has a water ice surface which is constantly 
bombarded by energetic ions. Atoms and molecules sublimated from the ice surface form 
a thin gravitationally bounded atmosphere. Suggestions that an oxygen atmosphere 
should accumulate above Europa’s water ice go back over 20 years [Wu et al. , 1978; 
Kumar and Hunten, 1982]. Atomic oxygen emissions from Europa was detected using 
the Hubble Space Telescope (HST) by Hall et al. [1995], who estimated the column 
abundances of molecular oxygen to be ~ 1.5 x 10 15 cm -2 and a surface pressure ~ 10 -a 
bar. These values are consistent with the equilibrium ice sublimation model, described 
by Johnson [1990]. 

Brown and Hill [1996] reported the discovery of sodium around Europa. Their 
conclusion is that this sodium originated from Europa’s surface. Sodium emissions are 
visible up to 25 Europa radii. Although sodium is a minor species, it should trace the 
distribution of dominant species. Thus, in addition to the exponentially decaying bound 
part there may be a large corona with an algebraic radial dependence. The best fit to 
the observations, performed by Brown and Hill [1996] favors ~ r -2 5 distribution of the 
neutral density (r is the distance from the center of Europa), which is somewhat less 
steep than the otherwise similar exponent (~ r -35 ) for Io sodium [Schneider et al ., 
1991]. A model for hot oxygen corona of Europa was suggested by Nagy et al. [1998]. 

Solar photoionization and particle impact ionization should produce an ionosphere 
around Europa. Galileo’s flybys on 19 December 1996 (E4) and on 20 (E6a) of February 
1997 were used for radio occultation and reconstruction of ionosphere profiles by Khore 
et al. [1997]. The typical distance between Galileo and Europa was about 1,500 km 
during the ingresses and about 4.000 km during the egresses. An additional distant radio 
occultation (E6b) took place on 25 of February 1997 when Galileo was about 2.8 xlO 6 
km from Europa. The latitudes of the occultations were close to —2°, —24° and —14° 
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respectively. These three occultations gave six ionospheric profiles which were averaged 
to give a plasma scale height of 240-440 km, with the maximum electron density of 10 4 
cm~ 3 near the surface. From these measurements Kl ore et al. [1997] inferred a neutral 
scale height of 120 km and a surface neutral density of 10 8 cm ~ 3 . Integrating this gives 
a column abundance of about 10 1 ’ cm~ 2 , which is consistent with the results of Hall 
et al. [1995]. However, there are very large uncertainties involved in the deduction of 
atmospheric densities from the ionospheric profile. It is based on just six occultations 
with significantly different latitudes and an impact ioiization rate which was most likely 
too low, just to mention two potential difficulties. One can see on figures in Kliore et al. 
[1997] that the electron density, obtained from the inversion of the radio-occultation 
data, often becomes negative. Thus, we will use the results of Kliore et al. [1997] as a 
starting guess, which will be modified, rather than a prescribed input for our model. 

The closest approach distance during Galileo’s tlyby of 19 December 1996 was 
688 km above the surface of the satellite. Successful magnetic field measurements 
were obtained during this flyby [ Kivelson et al., 1997]. Near the closest approach a 
perturbation of ~ 50 nT was observed, which was clearly not a random disturbance, but 
a coherent magnetic signature associated with Europa. Kivelson et al. [1997] explained 
the magnetic signature of Europa by assuming an intrinsic dipole field. Based on a 
magnetostatic vacuum superposition they have suggested an upper limit for Europa’s 
dipole which produces a maximum surface magnetic ield magnitude of 240 nT. In more 
recent papers Kivelson et al. [1998] and Khurana et al. [1998] suggested an induced 
dipole model as an explanation for at least part of the observed magnetic signature. 

During the E4 and E6a flybys Galileo plasma wave instrument observed strong 
whistler-mode noise. From the measured upper hybrid resonance frequency and the 
local electron plasma frequency Garnett et al. [1998 j inferred the electron density. 
They observed an electron density enhancement of cO to 100 cm -3 close to Europa. 
Unfortunately, the average charge state of the plasma ions is not known very well. If 
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we assume that it is about 2, then the results for the electron density of Gurnett et al. 
[1998] are close to the older model of Io plasma torus of Smyth and Combi [1988]. Thus, 
the Galileo measurements [ Gurnett et al. , 1998] imply a plasma density enhancement in 
the wake behind Europa 2-2.5 times the upstream values. 

While there is a wealth of models for ionospheres and magnetospheres of major 
planets and satellites, Europa until very recently has escaped the attention of the 
modelers. Prior to the Galileo mission very little was known about Europa, which 
resulted in only very crude models of the environment of the satellite. 

Wolff and Mendis [1983] considered the effects of the Jovian magnetosphere on 
the icy surface of Europa, Ganymede and Callisto; however they did not perform any 
simulation of plasma flow around the satellites. A decade later, Schreier et al. [1993] 
constructed a “zero-dimensional ' model of the Europa plasma torus. The focus of their 
work was the balance of the chemical species rather than plasma dynamics. Europa’s 
magnetospheric interaction was described by Ip [1996] from a physical point of view, 
while the main object of the paper was Europa’s exosphere. Some of the parameters 
used by Ip [1996], such as an acoustic Mach number for the upstream plasma of M = 2.7, 
seem to be inconsistent with our current understanding of Europa’s space environment. 

Recently, Saur et al. [1998] developed a more sophisticated model for Europa's 
magnetospheric interaction. However, the emphasis of their work was on the neutral 
atmosphere, for which they obtained some important results that agree with the 
observations of Hall et al. [1995]. At the same time, Saur et al. [1998] did not treat 
the magnetic field self-consistentlv, which diminishes the value of their results for the 
understanding of plasma flow around Europa and makes impossible an attempt to 
reproduce Galileo’s magnetometer data from their model. 

In a theoretical paper Neubauer [1998] considered general properties of sub-Alfvenic 
interaction of Galilean satellites with Jovian magnetosphere. He discussed the 
qualitative effects of ion pick-up, ion-neutral collisions and intrinsic magnetic field on 
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the plasma flow past a satellite and described several possible types of magnetospheric 
interaction. The quantitative understanding of the details of this interaction requires 
global three-dimensional MHD modeling. This is the main objective of our work; we also 
compare our model results with Galileo’s magnetometer and plasma wave instrument 
measurements. 


Model 


We describe the plasma environment of Europa with the single fluid ideal MHD 
equations. The conservative dimensionless form of these equations is 
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Here M is the local mass loading term (proportiona. to the neutral density), r\ is 
the ion-neutral collision momentum-energy transfer late per unit volume, \ is the 
recombination rate and po is the upstream plasma density. The derivation of these 
source terms from the Boltzmann equation can be foind in Gombosi [1991]. They are 
very similar to those used by Combi et al. [1998]. T le only difference is that in the 
present Europa model we included ion recombination. It is naturally treated as a second 
order reaction. Ion loss from recombination is accompanied by a production term whose 



role is to keep perfect balance for the upstream conditions. In our simulation we used 
Y = 8 x 10 -8 cm 3 s -1 based on Torr and Torr [1985]. This is the same value as used by 
Kliore et al. [1997]. This parameter was not adjusted through the simulation, but kept 
constant. 

Europa is approximated by a sphere 1500 km in radius. The boundary conditions 
on the surface were imposed bv utilizing cut cells [ DeZeeuw and Powell, 1993], which 
allow second order (piecewise linear) reproduction of the surface geometry. In our 
simulations we have studied a number of possible boundary conditions on Europa's 
surface. 

The simplest and the most numerically natural type of the boundary conditions is a 
perfectly conducting solid sphere. Under these assumptions, no plasma flow or magnetic 
flux are allowed to penetrate the body. However, this approach might not be physically 
appropriate. The conductivity of Europa’s ionosphere and surface ocean layer can be 
very high, but still finite, and thus, over millions of years of the planet existence, the 
external magnetic field would still diffuse into the planet interior. 

We can assume completely "absorbing” boundary conditions for Europa. The 
reason to allow the plasma to penetrate the surface of the satellite can be the following. 
Europa is surrounded by a neutral atmosphere, which is relatively dense near the 
surface. An impinging ion is likely to be lost in it, rather than diverted around the 
planet. Absorbing boundary conditions could also be explained by a recombination rate 
that becomes larger in the lower neutral atmosphere. A consistent model of ion-neutral 
interaction requires an ionosphere model for Europa to be coupled with the \1HD model 
of the magnetosphere. However, this kind of numerical simulation is beyond our present 
capabilities. The models of a superconducting sphere” and an “absorbing sphere” are 
clearly the two extremes bracketing reality. 

The system of MHD equations is solved using a higher order Godunov-type finite 
volume method on an unstructured Cartesian grid with adaptive mesh refinement. The 
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cell structure (“octree”) provides very good resolution in the interaction region, while 
allowing the total simulation box to be very large. Tl.e numerical aspects of the model 
and grid structure are described in DeZeeuw and Powell [1993], Powell [1994], Powell 
et al. [1995]. Convergence to the steady-state was accelerated by local time-stepping 
and optimal residue smoothing [can Leer et al., 1989] 

The reference frame we use throughout this article is a Cartesian system of reference 
in which the X axis is parallel to the corotation direction, Z axis is parallel to the Jovian 
spin axis and Y completes the right-handed system. 

The upstream parameters for Europa are not known exactly, so the input for our 
model was based on the Galileo magnetic field measurements [Kivelson et al., 1997] and 
on a model of the Io plasma torus [ Smyth and Combi 1988] derived from the Voyager 1 
plasma data. The following is a table of the parameters which we used: 

-plasma density 35 cm' 3 

-plasma temperature 90 eV 

-plasma speed 90.3 km/ s (corotational value) 

-mean molecular mass 24 amu 
-ratio of the specific heats 1.67 

-a uniform tilted Jovian magnetic field with components in the Cartesian frame of 
reference (63, -162, -422) nT (tilted 8° with respect to the upstream velocity). These 
parameters correspond to a sonic Mach number M = 0.7 and Alvfenic Mach number 
M a = 0.25. 

Results and discussion 

We have investigated a number of boundary corditions, mass-loading functions, 
ion-neutral charge exchange coefficients and intrinsic cipole fields trying to reproduce the 
Galileo measurements. In order to perform the comparison between the measurements 
and the model we have removed the global trend (c instant linear slope) from the 
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original magnetometer data [ Kivelson et ai, 1997], It is straightforward in principle 
to incorporate a constant slope (or even an arbitrary coordinate dependence) into the 
initial conditions for the magnetic field, however most of the large scale structure of the 
magnetic field is related to a time-dependent effect. This effect includes both Galileo 
orbital motion and Jupiter rotation around its spin axis which is 10° different from 
its dipole axis. Time-accurate calculation of the Galileo flyby is beyond our current 
computational capabilities and is not the purpose of the current paper. 

We show by the dashed line in figure 1, as a baseline, the results of vacuum 
superposition of a dipole field and Jupiter’s magnetic field, approximated by a constant 
vector. To get the best fit we have minimized the norm of the difference between the 
magnetostatic model and data. The resulting value for the dipole is (-0.1, 86., -108.) 
nT R E , which gives a magnetic field intensity of 270 nT at the surface at the magnetic 
pole. This dipole is within 15% of the value suggested by Kivelson et ai [1997]. 

The numerically simplest type of boundary condition is a perfectly conducting 
Europa. This type of boundary conditions was one of the two bracketing cases for Io 
used by Combi et ai [1998]. Figure 1 shows by dotted line the results of a simulation 
with this kind of boundary condition. In this case we did not use any mass-loading, 
but tried to fit the magnetic field measurements by changing the intrinsic dipole value. 
Our best results were obtained with a dipole moment of (16, 42, 57) R 3 E nT. This fit 
is no better than the magnetostatic superposition, but it has an important message 
- the dipole moment vector has almost reversed its direction as compared with the 
vacuum superposition. This tells us that the problem is very sensitive to the boundary- 
conditions. 

The major processes in Europa’s upper atmosphere leading to the production of 
new ions are photoionization and electron impact ionization of atmospheric neutrals. 
Kliore et al. [1997] estimated that the ionization frequency due to impacts is three 
times larger than that due to photons. This means that the ion production rates are 
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likely to be non-spherically symmetric and perhaps weighted toward the ram side. The 
actual angular dependence of the production rates is. however, unknown and cannot be 
inferred from available data. 

VVe started with the total mass-loading of 1.1 x 10 26 s _l distributed spherically 
symmetrically with an ion scale- height of 240 km, which are the parameters suggested by 
Kliore et al. [1997]. With mass- loading, ion-neutral charge exchange and recombination 
included we fail to obtain any reasonable fits to the data using a superconducting sphere 
for Europa, but had to switch to completely absorbing boundary conditions. Even with 
absorbing boundary conditions, the disturbance of the magnetic field in our model was 
too large and did not seem to have the correct character. Thus, we reduced the total 
mass-loading together with the scale-height. Our best fit - figure 1, continuous line - 
was obtained with a total mass-loading of 6 x 10 24 s~ l and an ion scale-height of 120 km. 
We assumed an ionization frequency of 1 x 10 6 s _l (same as Kliore et al. [1997]) and 
an ion-neutral charge exchange parameter of rj = 2.8 x 10 -8 cm 3 s -1 . These parameters 
are significantly different from those recommended by Kliore et al. [1997], but we point 
out the very large uncertainties that appear in their evaluation of the “best ionospheric 
profile”. A theoretical calculation Europa’s exosphere by Nagy et al. [1998] somewhat 
modified the earlier result of Kliore et al. [1997] and estimated the total escape rate 
of hot oxygen as 8.8 x 10 25 s _1 . Even this estimation is still significantly larger than 
our value for the total mass loading. In this simulation we used the dipole moment of 
(0,71,-64) R% nT which is approximately 40% smalle* than the value from the vacuum 
superposition. As one can expect, the mass-loading liminishes the required value of 
the intrinsic magnetic field. It does not seem feasible yet to set a lower limit different 
from zero on Europa’s magnetic field, because the details of the neutral atmosphere 
interaction with Jupiter's plasma torus are very uncertain and the boundary conditions 
have an important effect on the conclusions about Europa’s intrinsic magnetic field. It 
is quite possible that the mass-loading functions mor ; complicated than the spherically 
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symmetric ones used in this work can be responsible for the observed magnetic field 
signature. The future availability of plasma instrument data can significantly improve 
the precision of the model, since it will be possible to compare plasma densities and 
temperatures along the Galileo trajectory, which are most directly affected by the 
mass-loading parameters. The plots of plasma density, speed and temperature based on 
our model are shown in figure 2 and are in qualitative agreement with the measurements 
of Gumett et al. [1998]. 

The general structure of the plasma flow around Europa can be inferred from 
figures 3, 4 and 5. The density contours go from 40 cm -3 to 120 cm -3 with increments 
of 10 cm- 3 . Qualitatively, the density in the XY and XZ planes are similar, but one 
can see that in XZ plane the plasma tail behind Europa is much broader. This is a 
typical MHD effect explained by the anisotropy of the plasma flow with respect to the 
velocity direction. The magnetic field is mostly along the Z axis, which favors the flow 
around the obstacle in X 'Y plane and hinders the flow in XZ plane. The asymmetry 
between the two perpendicular slices of the three-dimensional solution is even more 
prominent in the velocity magnitude. The speed contours go from 88 km/s to 94 km/s 
in figure 4 and to 91 km/s in figure 5 in increments of 3 km/s. The lines with arrows 
are plasma stream-lines. They are not bent as much as they would be in the case of 
stronger obstacle, such as a conducting sphere. It is interesting to point out that the 
curvature of the stream-lines is different in the two cross-sectional planes due to the 
three-dimensional character of the flow. The speed increases faster and to a higher value 
along the Y axis than along the Z axis, which explains why density tail is broader in the 
Z direction than in the Y direction. 

Conclusions 

In this paper we presented the results of a three-dimensional MHD simulation of 
Europa magnetospheric interaction. Even in our relatively simple model we were able to 
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reproduce the main features of the Galileo magnetometer and plasma wave instrument 
measurements. Although we have used an intrinsic dipole moment for Europa in our 
model, it is possible that the satellite does not have one, and its magnetic signature can 
be explained without such a field. Many non-MHD effects, such as finite gyro-radius 
radius of ions, multi-fluid interaction, physical resistivity and viscosity, deviations from 
thermal equilibrium, cannot be self-consistently incorporated into our model yet. In our 
present model we have not included the effect of coronal-type mass-loading, which may 
be inferred from the measurements of Brown and Hill [1996] and the model calculations 
of Nagy et al. [1998]. It is clear, that this part of Europa’s atmosphere is very thin, 
and from our experience with the Io model [Combi et al ., 1998] the fine details of the 
mass-loading process are not overly important. In the future, however, the mass-loading 
far from Europa might play a critical role in interpreting the results of later Galileo 
flybys - [Kivelson, private communication ]. Study of these effects and comparison with 
coming Galileo data will be the basis of future work. 
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Figure Captions 
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Figure 1 . Magnetic field components: x - measured, dashed line - vacuum superposi- 
tion, dotted line - superconducting sphere with dipole and continuous line - absorbing 
sphere with mass-loading, ion-neutral charge exchange, recombination and dipole. 
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Figure 2. Density, speed and temperature profiles along the fly-by, for absorbing sphere 
with mass-loading, ion-neutral charge exchange, recombination and dipole. 
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Figure 3. Density; the contour lines are separated by 10 cm 3 increments. The lowest 
contour line (on the outside) is 40 cm' 3 . 
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Figure 4. Speed and streamlines in XY plane. The contour lines with 3 km/s incre- 
ments. The contours stretched in the Y direction correspond to 94 and 91 km/s , the 
outmost contour stretched in the X direction is 88 km/s. 
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Figure 5. Speed and streamlines in XZ plane. The contour lines with 3 km/ s increments. 
The “wings” correspond to 91 km/ s, the first closed :ontour is 8 8km/s. 
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Abstract. The development of the expanding atmosphere from the evaporating 
cometary nucleus has traditionally focused on observing and modeling the separate 
development of two distinct components, gas and dust, which are coupled dynam- 
ically with one another at distances out to a few tens of cometary radii. In the 
last decade or so, however, direct evidence from observations and suggestions from 
theory imply that the dusty-gas coma is a tightly coupled system where material 
is transferred between the solid and gaseous phase as an important integral part of 
the basic development of the coma. 

Comet Hale-Bopp (C/1995 Ol) was discovered far from the sun and is the largest 
and most productive comet, in the sense of release of gas and dust, in modern 
times. This has permitted observations to be made over an unprecedented range of 
heliocentric distance. This paper presents a review of a range of important issues 
regarding interrelations between dust and gas in comets, a description of the gas and 
dust environment for Hale-Bopp, and a summary of the preliminary results from 
Hale-Bopp which are relevant to these issues. Particular topics include dusty-gas 
models, dust fading and fragmentation, the role of dust and gas jets, the day/night 
distribution of gas and dust, and extended sources of gas. 

Key words: Comets, Hale-Bopp, Dynamics, Photochemistry 


1. Introduction 

Historically, the term comet has been given to the expanding atmo- 
sphere of neutral and ionized gas and dust surrounding a small icy 
parent body (the nucleus) which is in an eccentric orbit around the 
sun and which shows a head (or coma) and one or more tails. It is 
understood that evaporation of volatile ices from the comet's nucleus, 
the release of particulate material (dust), and the subsequent inter- 
action between dust and gas in its immediate vicinity (within a few 
nucleus radii) are responsible for the observed phenomena. For comet 
P/Halley we do have spacecraft images of the dust distribution in this 



2 


M R. COMBI ET AL. 


region. In general, however, the innermost coma server as an unresolved 
black box source for the gas and dust observed remotely outside this 
region. The detailed connection between the icy parent body and the 
black-box source of gas and dust remains clouded. There is still an 
incomplete picture regarding the overall nature of the conditions and 
processes occurring in the innermost coma. As will be discussed in detail 
below, recent evidence has blurred the traditional separation between 
the dust and gas components of the coma, indicating that not only are 
they dynamically coupled to one another but possibly together form a 
tightly coupled system. There is strong observational evidence for the 
production of extended non-nuclear sources of gas species which might 
come from dust or icy grains, and there are theoretical suggestions of 
recondensation of gases back into icy grains. 

Comet Hale-Bopp (C/1995 01) was discovered far from the sun and 
is the largest and most productive comet, in the sense of release of gas 
and dust, in modern times. This, combined with continued advances 
in ground-based and near- Earth-based instrumentation, has permitted 
an unprecedented range and quality of remote observations of gas and 
dust to be obtained. Because the geocentric distance for Hale-Bopp 
was never very small, this inner region is typically unresolved, and 
therefore continued progress in our understanding of the important 
innermost coma requires a combination of observational evidence of gas 
and dust exiting this region and theoretical modeling which accounts 
for the basic physics and chemistry. Observations of :omet Hale-Bopp 
which provide new insight into the interrelations betvreen gas and dust 
in the coma are summarized. Theoretical modeling developments in 
the study of dusty-gas dynamics axe also discussed. Particular top- 
ics include dusty-gas models, dust fading and fragmentation, the role 
of dust and gas jets, the day/night distribution of gas and dust, and 
extended sources of gas. 


2. The Gas-Dust Environment of Comet Hale-Bopp 

Most coma modeling has been performed either for rather faint, to 
moderate, to so-called bright (P/Halley) comets. The very large pro- 
duction rate of Hale-Bopp introduces some interesting physical effects 
which axe usually not very important in other comets and are often 
totally ignored in the interpretation of data. The first obvious differ- 
ence is the large collisional coma. A typical definition for the collisional 
zone of a spherical coma is that distance where th* local collisional 
mean free path is equal to the distance from the cenler of the nucleus. 
For example Whipple and Huebner (1976) give 
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For the conditions of a Halley-class comet at 1 AU, a total gas pro- 
duction rate, Q ga s — 7x 10 29 s~ [ , a nominal collision cross section of a = 
3 x 10“ 15 cm 2 , and an outflow velocity of about 0.9 km s _1 , the collision 
zone radius is about 1.9 x 10 4 km. For comet Hale-Bopp near perihelion 
the total gas production rate is close to 10 3l s _1 . All else being equal, 
this would yield a collision sphere of about 2.7 x 10 5 km. The primary 
heating agent of the neutral coma is the collisional thermalization of 
the energy imparted to hot H atoms produced in the main branch of 
water photodissociation. In a Halley-class comet, because of the 20-to- 
1 mass ratio between H atoms and the average heavy molecule, the H 
atoms actually begin to decouple collisionally at 10 % of the typical 
collision zone radius, producing a moderating effect on photochemical 
heating which is dependent on the magnitude of the gas production 
rate. This was first treated in the context of hydrodynamic models 
using heuristic approximations (Ip 1983; Kitamura 1986), escape prob- 
abilities (Huebner and Ready 1983), and multifluids for pre-thermal 
and thermalized H (Marconi and Mendis 1984). Later iterative Monte 
Carlo-type kinetic models for H atoms were developed to explicitly cal- 
culation the collisional efficiency (Combi 1987; Bockelee-Morvan and 
Crovisier 1987; Combi and Smyth 1988; Ip 1989). The combination of 
the variations of photochemical heating rate and the overall gas pro- 
duction rate with heliocentric distance is responsible for the variation 
of the terminal outflow speed of molecules and radicals in the outer 
coma. The large collision sphere implies that comet Hale-Bopp should 
have had a much larger region of photochemical heating resulting in 
higher gas kinetic temperatures and larger outflow speeds. A second 
related effect of the large gas densities is that infrared radiative cooling 
is effectively shut-off for the entire collisional coma (Bockelee-Morvan 
and Crovisier 1987), thereby further enhancing the effect of the photo- 
chemical heating on the coma temperature and outflow speed. 

The last important effect of the large coma densities on the expected 
behavior of the coma is the non-negligible opacity of ionizing and dis- 
sociating solar ultraviolet radiation into the inner coma. Marconi and 
Mendis (1984) had calculated a ID multifluid spherical coma model 
for comet Halley conditions, with a self-consistent calculation of the 
effect of UV opacity in three UV wavelength bands which are primarily 
responsible for the main photodecay branches of water: 
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1: A < 1860 A: H + OH 
2: A < 1450 A: H 2 + 0( 1 D) 

3: A < 984 A: ionization (mostly H 2 0 + ). 

The half-power points of their radial profiles cf UV attenuation 
serves a good indicator of the typical distance for shielding the inner 
coma from incoming UV solar radiation in each wavelength band. If 
their results are scaled from the production rate in Halley to that of 
Hale-Bopp at 1 AU, shielding distances of 600, 1100, and 3300 km in 
each of the three wavelength bands, respectively, are found. This means 
that in the case of a spherical coma for comet Hale-Bopp ionization will 
be effectively shut off within ~3300 km of the nucleus (and in a cylindri- 
cal shadow extending to ~60,000 km antisunward given the size of the 
Sun l s disk). Photodissociations will similarly be shut-off within about 
800 km of the nucleus, and the ratios between the various branches will 
be modified within distances of ~1000 km. For a non-spherical coma, 
for instance where the gas production is skewed primarily sunward 
and in narrow jets, local effects could be even stronger because of the 
higher concentrations of gas density. The increased opacity could have 
consequences on extracting parent production rates or interpreting spa- 
tial profiles obtained from very high spatial resolution observations of 
daughter radicals near the nucleus. 

Using a ID-spherical hybrid kinetic/dusty-gas-hydrodynamic calcu- 
lation with a fairly simple composition, the gas velocity and tempera- 
ture have been calculated for a range of heliocentric distances for condi- 
tions appropriate for comet Hale-Bopp. This formulation accounts for a 
water-driven photochemistry, which should be dominant even for fairly 
large heliocentric distances. The kinetic portion of the calculation is a 
Monte Carlo model that explicitly accounts for the collisional efficiency 
of hot H atoms from the main water dissociation, and is solved iterative- 
ly with the ID-spherical steady-state equations for mass, momentum, 
and energy assuming an ideal gas. The detailed parameters (lifetimes, 
branching ratios, collision cross sections) have been discussed at length 
already (Combi and Smyth 1988; Combi 1989; Ccmbi and Feldman 
1993). For these models a nucleus radius of 25 km is adopted and the 
gas coma is assumed to be composed of 80% H 2 0, If .5% CO, and 3.5% 
C0 2 . The CO and C0 2 are included mainly to raise the mean molecu- 
lar mass to the correct level (about 20 amu). All of the dust is assumed 
to have a mean size of 0.65 mm. Past calculations fo' P/Halley (Combi 
1989) produced reasonable results for both the velo :ity of the gas and 
the dust particles using only a single dust size of 0.65 /im and a density 
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Table I. Hybrid Kinetic Dusty 
Gas Hydrodynamic Model. 

r Q(H 2 0) 

(AU) ( 10 28 s” l ) 



40. 

2.5 

70. 

2 

140. 

1.5 

200. 

1 

800. 


of 1 g cm -3 which gave equivalent results to using an entire more real- 
istic distribution of dust sizes and size dependent particle bulk density 
(Gombosi et ah 1986). For the calculations presented here a dust-to-gas 
ratio of 5 times that of Halley was adopted, as is indicated in photo- 
metric measurements of Afp by Schleicher et al. (1997). As will be 
discussed below, Afp is only a good indicator of the abundance of vis- 
ible wavelength (half-mic ron) sized particles, and is only related to the 
real dust-to-gas ratio given some the particle size and terminal velocity 
distributions. The calculation accounts for the UV optical depth of the 
water photodissociation rate which is effectively shut off in the inner 
coma for distances ranging from 800 km at 1 AU to only 35 km at 
3 AU. Table I gives the assumed parameters used for the particular 
calculations. 

Figure 1 shows the variation of the gas radial outflow velocity and 
the gas kinetic temperature as a function of distance from the center 
of the nucleus for the model calculations. The temperature shows the 
typical drop to very low values at a few hundred km caused by the 
dominance of the pressure-driven expansion. The precise form of the 
drop is a balance between the expansion-cooling and dust-gas heating in 
the inner coma and photochemical heating in the outer coma. Because 
of the large production rate, rather larger outflow speeds for parent 
molecules are expected to be found for Hale-Bopp than for Halley at 
comparable heliocentric distances. These are produced by increased 
photochemical heating efficiency resulting from the more complete ther- 
malization of hot H atoms, as well as by increased gas heating by the 
dust in the inner coma caused by the large dust-to-gas ratio. 

These model results compare favorably to early observations from 
which outflow velocity and temperature of the coma have been deter- 
mined. Biver et al. (1997. 1998) derived mean coma outflow speeds and 
rotational temperatures from their mm and sub-mm observations of CO 
and CH 3 OH at IRAM and .JCMT corresponding to beam aperture full- 
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Figure 1. ID-Spherical Hybrid Kinetic/ Dusty- Gas Hydrodynamics. The results of 
ID-spherical dusty-gas hydrodynamic model for comet Hale-B<>pp where the photo- 
chemical heating rate is calculated by a kinetic Monte Carlo s mulation. 


width-half-maxima (FWHM) in the range of 10 to 26 arcsec. At least in 
the inner coma, which is primarily sampled by these measurements, the 
large collision rates discussed above, should maintain a reasonable equi- 
librium between the gas kinetic and rotational temperatures. We have 
integrated the mean radial outflow velocity and mea i gas kinetic tem- 
perature from the results obtained from the 1-D mod d over a Gaussian 
beam of 14 arcsec (the intermediate size for the observations) at the 
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appropriate pre- and post- perihelion geocentric distances and fitted 
them to power-law variations with heliocentric distance for comparison 
with the results of Biver et al. They found pre- and post-perihelion 
velocities in km s -1 of 1.18 r -0 - 44 and 1.07 r“°' 37 , and rotational tem- 
peratures of 116 r _1 24 and 93 r _a83 in Kelvin. Here, r is the heliocentric 
distance in AU. For the corresponding model values we find pre- and 
post- perihelion velocities of 1.01 r" 032 and 1.02 r~ 034 and temper- 
atures of 114 r“ 1,6 and 117 r~ 16 . The small differences between pre- 
and post- perihelion results from the model are due to the differences in 
geocentric distance for corresponding observations before and after per- 
ihelion. The agreement between model and observations for the mean 
outflow speed is quite good. Perhaps this should not be too surpris- 
ing because the same calculation was successful in explaining a num- 
ber Doppler resolved observations for comets P/Halley and Hyakutake 
1996B2 (Combi 1989; Smyth et al. 1995; Combi and Cochran 1997). 
Despite the fact that it is not clear that the mean rotational tempera- 
ture extracted from the observed aperture-integrated spectrum should 
necessarily yield the same temperature as the mean value of the gas 
kinetic temperature integrated over the same aperture, the agreement 
between temperatures and their variation between the observations and 
the simple model calculation is nonetheless quite good. 

An independent confirmation of the model outflow speeds at large 
distances is found in a comparison with the analysis of the OH line 
profile observed with the Nangay radio telescope (HPBW of 3.5 by 19 
arc min) by Colom et al. (1998). Using a model which assumes an OH 
velocity of 1.05 km s -1 , water and OH lifetimes at 1 AU of 8.2 xlO 4 
s and 1.1-1. 6 xlO 5 s (the range coming from the heliocentric velocity 
dependence) they find mean water outflow velocities of 0.7 km s -1 for 
the comet at 3 AU and 2.2 km s _l at 0.93 AU. This is in excellent 
agreement with the outflow speeds at large distances from the comet 
at both heliocentric distances in the model calculation. 

Measurements of gas species are typically interpreted using a model 
which adopts (by necessity) an outflow speed. The radial outflow speed 
for Hale-Bopp clearly varies with both distance from the nucleus and 
overall with heliocentric distance. Based on our calculations both vari- 
ations can be expected to be much more pronounced than had been 
found for Halley (Combi 1989) because of the enhanced photochemical 
heating efficiency in Hale-Bopp. Based on the comparison above with 
velocities and temperatures extracted from observations, we believe 
that the model velocities and temperatures shown here are reasonably 
representative of the conditions in comet Hale-Bopp. The velocities 
might be especially useful in the analysis of observations where the 
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dynamical conditions in the coma are needed. These are given in Table 
II. 

Observations of the coma of comet Hale-Bopp indicate that the ejec- 
tion of dust into the coma is primarily, perhaps exclusively in the sun- 
ward direction. (Rauer et al. 1997; Hayward et al. 1997; Weaver et al. 

1997) . Near perihelion a well-developed system of dust haloes on the 
sunward side of the nucleus indicates* that dust is ejected into a nar- 
row jet from a primary active area when that area is illuminated by 
sunlight (Lederer et al. 1998; Larson et al. 1998). Gas seems also to be 
primarily produced from this active area, although it appears to remain 
active at a reduced level even when rotated into the apparent night side 
of the nucleus (Samarasinha et al. 1998; Larson et ai. 1998). It is clear 
then that, although ID spherical models are valuable for determining 
the nature of global average gas and dust dynamics, multidimension- 
al models are necessary to understand many important attributes of 
comet Hale-Bopp. 

We present here the results of a three-dimensional calculation for 
the dusty- gas dynamic flow in the vicinity of the nucleus for an active 
area superimposed on a uniform background. The calculation is per- 
formed on an unstructured multiscale Cartesian octree grid using adap- 
tive mesh refinement. This octree structure is a hierarchical cell struc- 
ture, based on multiple generations of cubic parent :ells which can be 
split into 8 child cells (see DeZeeuw and Powell 19)2). The adaptive 
refinement of the mesh follows the gradients of the dynamical variables 
of the calculation and naturally results in the smallest cubes in regions 
where the flow gradients are the highest. Such a scheme makes optimal 
use of computer resources allowing problems with disparate scales and 
complex geometries to be modeled globally (i.e., a large simulation vol- 
ume), while still retaining sufficient resolution for following small fea- 
tures locally , such as capturing shocks. The general numerical method 
has been applied to the MHD equations for the comet solar- wind inter- 
action (see Gombosi et al. 1996; Haber li et al. 199" 'a&b; Gombosi et 
al. 1998), and now to neutral dusty-gas flows without B-fields (Kabin 
et al. 1998). 

Our method differs from a series of limited three-dimensional dusty- 
gas calculations by Crifo and coworkers (Crifo et a.. 1995, 1997, and 

1998) in fundamentally important ways. With then model they have 
already investigated a number of complicated emission conditions and 
geometries, including, weak and strong comets, miltiple jets, many 
dust sizes, and non- spherical nuclei. Their calculetion is performed 
on variable- sized, but fixed, spherical-polar grid. Because of the natu- 
ral singularity in the spherical-polar grid at the pol?, their solution is 
limited to regions away from the pole. Our method lses nested Carte- 
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Table II. ID Spherical Hybrid Kinetic/Hydrodvnamic 
Gas Outflow Velocities. 
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d = distance from the center of the nucleus in km. 
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sian cells, avoiding the geometrical singularity limitation and permit- 
ting fully global modeling adapting the cell structure to the calculation 
automatically, and with no spatial limitations. 

Briefly the calculation solves the standard dusty-gas Euler equations 
for mass density, momentum and energy for a single-fluid gas and the 
equations for mass density and momentum of a number of dust particle 
sizes. For a detailed discussion see the review by Gombosi et al. (1986). 
The coupled system can be written in Cartesian coordinates as: 

Op _ . . 5p 

^ +v( " u)= « 


+ p { u • V)u + Vp = -F 
at 


1 dp + — L_ (u .V)p+^— p(Vu) = -Q 


7 — 1 dt 7-1 
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dpi , „ , . & Pi 

sr + v (ftVt) = It 
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f)"V 

Pi-^+P, (v , V )v , = F < i = 1... V 

where p is the gas mass density, u is the gas velocity, p is the gas 
pressure, vi and pi are the velocity and mass density for dust particles 
of radius a,. The terms on the right hand sides of the equations are the 
various source terms. & is the gas production source. ^ is the dust 
production source. F is the gas-dust drag force. Q is the dust-gas heat 
exchange rate. Then we can write the source terms <s 
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Here C v is the gas heat capacity at constant pressure and the rest 
of the coefficients can be defined under the assumption of diffusive 
reflection as: 



The dust bulk density, pd , is the mass density of a particular dust 
particle. The gas temperature is T. The dust temperature, T*, is con- 
stant for a particular dust particle size in equilibrium with solar radi- 
ation. C p is the gas specific heat at constant pressure and the various 
other quantities are given assuming the standard free molecular condi- 
tions for gas-dust (Probstein 1968 with minor corrections by Kitamu- 
ra 1986). These are s z , the gas-dust relative Mach number, C d\ the 
dust-gas drag coefficient, T[ ec , the recovery temperature, R/, the heat 
transfer function, and St/ the Stanton number. The equations are pre- 
sented here in terms of the primitive variables (density, velocity, and 
pressure); they are rewritten in terms of the conserved variables (mass, 
momentum, and energy) for the purpose of solving them numerically. 

We first present here a calculation for a single concentrated gas/dust 
jet with uniform production over the rest of the nucleus. The nucleus 
is taken to be a sphere with radius of 28 km. Based on rough estimates 
relevant for comet Hale-Bopp at 1 AU, we have set 30% of the gas 
and dust production to come from a single active area with a Gaussian 
in angle having a width of 10° (half-angle), and the remaining 70% 
spread uniformly over the nucleus. The total gas production rate of 
10 31 s -1 with a mean molecular mass of 20 amu is adopted. Six sizes 
of dust particles with radii of 0.1, 1, 10, 100, 1000, and 10000 gm 
(the latter three corresponding to 0.1, 1 and 10 mm) are adopted with 
a dust- to-gas mass ratio of 5 (Schleicher et al. 1997) and a number 
distribution varying as a -3 where a is the particle radius. The precise 
value of the dust-to-gas mass ratios depends critically on the extension 
of the population to very large particles, which contain most of the 
mass, and the functional form of the size distribution which is not 
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well-specified. For example, essentially the same distribution of small 
particles for comet Halley yielded total dust-to-gas mass ratios varying 
from 0.3 (Gombosi 1986) to greater than 1 (Crifo 1991). The actual 
difference was only in the extension of the distribution to the very large 
particles (cm-sized) by Crifo, which effect neither the gas dynamics, nor 
contribute to ground-based observations of dust continuum. 

For the model calculations shown here the entire simulation box 
extended to ± 50 times the nucleus radius (Rn), or to 1400 km from 
the nucleus. A total of 220,000 cells were included in the simulation 
and the cells ranged in size from 0.78 to 10 Rn. The adaptive mesh 
structure is shown for the full simulation box and for i magnified region 
near the nucleus in Figure 2. As in the models of Crifo and co-workers 
each dust size is a separate fluid, however in our case the full multi- 
fluid calculation is carried out in the whole simulation volume. The 
full multifluid calculations of Crifo and co-workers which include many 
more dust sizes (up to 44 logarithmically spaced s zes) are typically 
carried out in a restricted volume near the nucleus, and are extend- 
ed to large distances outside the subsonic region us ng a steady-state 
expansion calculation. In our full calculation each dust size introduces 
four additional equations and unknowns to be solved in every cell. 

There are a number of important and interesting results from our 
dusty-gas calculation. Plate la-b (color section) shows color contour 
plots of the number density of the gas with superimposed velocity 
streamlines and the gas speed. Note that the gas flow lines originating 
from the active area expands tangentially near the nucleus but that by 
a several R^ the flow is essentially radial and the £ symmetry due to 
the active area is frozen-in. Plates 2a-b and 3a-b (color section) show 
similar pairs of plots for two of the sizes of dust with radii of 1 fj,m 
and 1 mm. In all cases the color table for the speed plots is the same. 
Colors for the density are relative because of the large differences in 
dust bulk mass density from size to size and compared with the gas. 
We selected these two particular dust sizes because both the density 
and speed variations are typical for the nearby dust sizes. The plots 
for the 1 fim particles are very similar to those for the 0.1 and 10 fi m 
particles. Those for the 1 mm particles are similar to those for the 0.1 
and 10 mm particles, although the speeds of course are somewhat larg- 
er and smaller, respectively, for the two groups. For all dust sizes the 
dust jet remains more concentrated over the active area and the gas 
squirts laterally to the side of the dust jet, making tl e gas distribution 
more uniform over the day-side hemisphere than ttu dust. This is not 
unexpected. 

The dust velocities for the small dust sizes (0.1, 1, and 10 mm) follow 
the gas flow quite closely. Because of the dust mass loiding, the gas and 
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Figure 2. Two Views of Adaptive Mesh for Dust-Gas Jet Calculation. Above in (a) 
the whole simulation volume is shown extending to 50 nucleus radii (1400 km) from 
the nucleus. Below in (b) the cell structure around the nucleus is shown, where the 
nucleus corresponds to the black circle at the center. 


small-dust speeds are slowest in the middle of the dust jet and fastest 
throughout the uniform background. Owing to the large gas densities 
the smaller dust sizes follow the gas flow very closely everywhere, and 
only the larger dust sizes show the conical type structure which has 
been seen even in the early 2D- axisymmetric dusty-gas jet modeling 
(Kitamura 1986; Korosmezey and Gombosi 1990). Furthermore, the 
larger dust sizes achieve the highest speeds in the jet (i.e., to the left in 
the Plate 3b), where the gas densities are the highest, and somewhat 
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lower speeds away from the jet (i.e. to the left in Plate 3b), the opposite 
behavior of the small particles. 

The most important implications from the model calculations result 
from the overall large dust speeds of all particles, compared with past 
calculations of moderately bright, Halley-class, comets. All of the parti- 
cles in the optical size range, i.e. those that make the dominant contri- 
bution to optical CCD imaging of the dust coma, and of course smaller 
particles, are accelerated to ~90% or more of the gas velocity. There is 
only a small range in speeds. Even the 10 /im particles are accelerated 
to ~75% of the gas speed. The uniform speeds of particles in the optical 
size range may be responsible for the sharpness of the observed dust jets 
in Hale-Bopp and for the rough spatial coincidence (or at least roughly 
similar spiral pitch angles) of the spiral dust jets and gas jets in the 
CCD-images (Larson et al. 1998; Lederer et al. 1998; Samarasinha et 
al. 1998). In Halley, the small grains suggested to be the source of the 
gas jets (sub-^zm) were accelerated to nearly the gas velocity, but the 
larger optical-size grains (^m sized) were not. The combination of dif- 
ferent outflow speeds with nucleus rotation produced a velocity sorting 
effect between gas (CN) and dust jets (Samarasinha, private commu- 
nication). In Hale- Bopp near perihelion, our calculation shows that 
all the optical and smaller particles are accelerated to nearly the gas 
speed. Therefore, there should be little velocity-sorting between optical 
and smaller particles, resulting in sharp dust jets, and gas jets following 
the same spiral arcs as the dust jets. We have also run a model with 
85% of the gas production and all of the dust production coming from 
one active area, and the resulting large dust velocities for optical- sized 
particles are similar (see below). 

Only the largest particles show the typical fall-oif in velocity with 
increasing particle size, however, the velocities are still much larger 
than the usual published values that are based on comets with much 
lower gas production rates (Crifo 1991; Gombosi et al. 1986). This may 
have important consequences on the estimate of tota 1 dust-to-gas mass 
ratios extracted from long- wavelength observations v r hich are sensitive 
to detecting most of the dust mass that is present in nm-sized particles 
(Senay et al. 1998). In Table III we show values of the gas and dust 
particle velocities at different angular locations compared with the jet 
near the outer portion of the simulation box at a distance of 1300 km 
from the center of the nucleus. 

Only preliminary reduction and analyses have b?en performed on 
the data needed to properly characterize the realistic balance between 
the separate components of gas and dust in an active jet and in a 
uniform background. The two cases we have examined, one with 30% 
of the production in the jet and one with 85% of the production in the 
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Table III. ID Gas and Dust Terminal 
Speeds in a 3D Dust-Gas Jet 


Terminal Speeds in m s 


Species 

In Jet 

45° 

o 

0 

135°- 180‘ 

Gas 

718 

734 

748 

810 

0.1 ^ni 

698 

718 

727 

795 

1.0 /xm 

647 

668 

669 

741 

10 

549 

556 

555 

603 

0.1 rmn 

394 

370 

370 

355 

1.0 mm 

212 

182 

182 

148 

10 mm 

84 

72 

73 

48 


jet, yield very similar large dust particle speeds. Therefore, the large 
dust speeds can be used in the next iteration of the process of trying 
to understand the temporal and geometric aspects of rotating dust and 
gas jets in the coma from the observations. They should also be used in 
the interpretation of long wavelength observations which detect mainly 
the large mm-sized particles which are the dominant contributors to 
the global dust mass production rate. 

For comparison with the first model results shown in Plates 1-3, 
Plate 4a-b (color section) shows the gas number density and dust bulk 
mass density for the 1-^m optical- sized particles in the vicinity of the 
nucleus in the case of a strong jet where 85% of the gas production 
and 100% of the dust production is in the jet and only 15% of the gas 
production is distributed uniformly around the nucleus. For this more 
confined distribution of material at the surface of the nucleus the dust 
jet remains more focused but the gas distribution broadens, expanding 
to 2 to 3 times its original Gaussian width. Some evidence from the 
observed images indicates that the observed jets may be narrower than 
even this model where 85% of the production is in the jet. This could 
indicate a role for some kind of collimation due to the spatial topology 
of the surface, which is not in our current model. Such collimation could 
be caused by the active area that produces the jet being recessed into 
the surface as in a crevasse or a crater. The effects of surface topology 
and local solar insolation which has possible influence on the strength 
and appearance of gas and dust jets have begun to be addressed by 
Crifoet al. (1998). 

Boice et al. (1998) have also recently presented results of a ID- 
spherical dusty coma model of comet Hale-Bopp at several heliocen- 
tric distances which includes the surface layer as part of the calcu- 
lation. Their calculation accounts for gas emission from the surface, 
from transport through a porous surface layer, and from release of gas 
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from dust in the coma. Separate fluids are included for electrons, heavy 
neutral gas, fast hydrogen, and dust. 


3. Evidence from Observations 

3.1. Extended Sources of Gas 

The discovery of an extended source of formaldehyde (H 2 CO) in comet 
P/Halley by Eberhardt et al. (1987) using data from the Neutral Gas 
Mass spectrometer on the Giotto spacecraft provided conclusive evi- 
dence of a substantial source of a gas species from unseen parent grains. 
There was not a sufficient amount of another heavier parent molecule 
to account for H 2 CO which was present at the level compared 

with water. That H 2 CO is produced by an extended source in comet 
Hale-Bopp has also been suggested by Womack et al. (1997) and Wink 
et al. (1998). Radio observations of HNC by Biver et al. (1997) and mm- 
wavelength observations of CH 3 OH by Womack et al. (1997) indicate 
that these species are also produced by an extended source. However, 
whether they are produced by grains or molecules (e.g., by photodisso- 
ciation) is not yet known. In addition, the H 2 CO itself was also found 
to serve as an important extended source for a subs* antial fraction of 
cometary CO (Meier et al. 1993; Samarasinha and Eelton 1993). This 
has since been seen in infrared observations of a numoer of subsequent 
comets including most recently Hyakutake and Hale-Bopp (DiSanti et 
al. 1998). From radio observations of the OH radical in Hale-Bopp, 
Colom et al. (1998) report evidence for water sublimation from grains. 

3.2. Dust and Gas Jets, Haloes and Shells 

One of the most interesting aspects of the early analysis of observations 
of Hale-Bopp is the dominance of gas jet activity. Shor :ly after discovery 
at nearly 7 AU it was found that CO, probably the dominant species at 
that time, was being produced in a naxrow, well-focised, day-side jet 
(Jewitt et al. 1996). Subsequently dust jets in various configurations 
were seen to dominate continuum images (Kidger et al. 1996; Weaver 
et al. 1997). Near perihelion in addition to the wo 11- developed day- 
side dust haloes, which were apparent to anyone loo dng at the comet 
through even an amateur telescope, narrow-band CCD imaging showed 
a set gas jets/haloes. The dust sources were only substantially active 
on the day side of the nucleus. This created a pattern made by one 
or more spiral jets which were not visible on the ni^ht side. The gas 
jets, however, continued activity on the night side, bu: at reduced level. 
Therefore, the gas jet spiral pattern was generally more continuous. 
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The CN radical, along with C 2 , was one of the earliest spectroscop- 
ic detections made in comets (Huggins 1882). Coming from the pre- 
sumption that CN was produced by the photodissociation of some sta- 
ble parent molecule, early studies (Potter and DelDuca 1964; Jackson 
1976) lead to the suggestion of various candidate species (HCN, C 2 N 2 , 
CH 3 CN, HC 3 N, etc.). The study of a large sample of both radial spatial 
profiles of CN (Combi and Delsemme 1980) and sunward and antisun- 
ward profiles distorted by radiation pressure acceleration (Combi 1980) 
lead them to conclude that the observational evidence was consistent 
with the production of CN by the photodissociation of HCN which had 
been discovered by radio emission (Huebner et al. 1974). However, their 
results were based on modeling the parent molecule outflow speed using 
the a heliocentric distance dependent velocity law of 0.58 r ~ 0 ' 6 km s -1 
(Delsemme 1982) which is now known to be too small. Later compar- 
isons of production rates from radio observations of HCN and optical 
observations of CN (Bockelee-Morvan and Crovisier 1985) showed that, 
although for some comets the production rates were roughly equivalent, 
for others the CN production rate was too large. This indicated mul- 
tiple sources for CN with different relative compositions, varying from 
comet to comet. 

The discovery of jets of CN in comet P/Halley by A’Hearn et al. 
(1986) lead them to the conclusion that the CN seen in the jets (at 
least) was produced from very small CHON (organic particulates rich in 
compounds containing carbon, hydrogen, oxygen, and nitrogen) grains. 
They reasoned that only a grain source would remain tightly entrained 
in the collimated flow. Although it was clear that most gas species and 
dust were produced from the same major active areas on the surface of 
the rotating nucleus (Millis and Schleicher 1986), the non-coincidence 
of CN and dust jets was explained as a velocity sorting according to 
grain size. The dust continuum was dominated by micron-sized grains 
which are accelerated only to some fraction (0.5 or less) of the gas 
velocity (Gombosi et al. 1986), whereas the CN in the jets was produced 
by very small CHON grains, which are accelerated to nearly the gas 
speed. Combi (1987), on the other hand, suggested that a localized 
source of parent molecules might remain collimated enough without 
the requirement to invoke grains. Parent molecules obviously travel at 
the gas speed so the non-spatial coincidence would be same whether 
produced from molecules or sub-micron CHON grains. A subsequent 
analysis of the CN-jet images by Klavetter and A’Hearn (1994) showed 
different radial distributions for CN within the jet compared to the 
CN between the jets, pointing again toward a CN-jet source (grains) 
with a length scale of ^8000 km, similar to that observed for H 2 CO. In 
addition to CN and C 2 jets (Schulz and A'Hearn 1995), Cosmovici et 
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al. (1988) also reported jets in images of C3, as well as in a combination 
image of [OI] plus NH2 emissions in comet P/Halley. 

For comet Hale-Bopp pre-perihelion observations of HCN (Biver et 
al. 1997) and CN (Schleicher et al. 1997 and Rauer et al. 1997) indi- 
cate there is a rough equivalence of the HCN to CN production rate 
for heliocentric distances larger than about 4 AU. By 3 AU the relative 
CN production rate increases over HCN by about 60%, and although 
the published data are incomplete at the time of zhis writing, indi- 
cations are that the ratio is approaching 2 near perihelion. Such a 
behavior is consistent with the gradual activation of a second source of 
CN at smaller heliocentric distances. Rauer et al. (1997) suggest possi- 
ble sources in grains which might be activated by 3 AU, or additional 
parent molecules (CH3CN and HCN). 

Fitzsimmons and Cartwright (1996) and Wagner and Schleicher 
(1997) suggest that the appearance of CN spectroscopically at greater 
than 6 AU from the sun points to the same mechanism for at least a 
large fraction of the CN seen at small heliocentric distances. Therefore, 
a reasonably consistent picture emerges for a combination of sources. 
The major one is probably the photodissociation of some set of more 
volatile parent molecule species including HCN which are active at large 
heliocentric distances. Then at smaller distances (3 AU) an addition- 
al source becomes activated which could be either less volatile parent 
molecules, small vaporizing grains, or some combination. The variation 
of the ratios of the various sources from comet to comet, with time 
(heliocentric distance), and possibly with location on the nucleus, could 
explain many of the difficulties in tying together a simpler consistent 
picture with HCN alone (Bockelee-Morvan and Crovisier 1985). 

Regarding the identity of other photodissociation parents for CN, 
Festou et al. (1998) are suggesting C2N2, which has been already on 
candidate lists for a long time. They find that the a combination of 
a lifetime of 3.5 x 10 4 seconds at 1 AU and a CN velocity of 1-2 km 
s^ 1 can explain the observed characteristics of production scale length 
and the Doppler-broadened velocity profiles. However, caution must be 
taken that model analysis must account properly for the expected large 
parent molecule outflow velocities (2 km s _1 at 10 5 km) and collisional 
quenching of the ejection velocity received by daugiter radicals upon 
dissociation of their parents. These attributes of comet Hale- Bopp have 
been discussed above. 

Lederer et al. (1998) have reported gas jets in all of the major visible 
neutral species, CN, C 2 , C 3 , NH, and (most impor;antly) OH during 
the two-month period surrounding perihelion. Wink et al. (1998) have 
presented Doppler- resolved radio observations during the same period 
of CO, HCN, CS, and H 2 CO, which indicate produ:tion these species 
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into localized jets rotating with the nucleus. The CO jet is active on 
both the day and the night side and appears to be located at the equa- 
tor. The average line-of-sight velocity of CO varies periodically with the 
11. 3- hour rotation period of the nucleus and has an amplitude of about 
0.2 km s _l . This is not necessarily indicative of the ejection speed of 
gas from the jet, but the integrated deviation from a stationary mean 
composed of the coma average* with the velocity dispersion of the jet 
projected along the line of sight. Wink et al. find that the HCN emission 
comes from the same location as CO but is weakened in the night side. 
The CS emission changes in a similar manner but appears to come from 
an active area in a different location in the northern hemisphere. The 
H 2 CO emission appears in phase with CS but is systematically more 
blue-shifted. Because FOCO is not a parent species (likely produced by 
grains), the extra velocity shift may not be surprising. 

Larson and Hergenrother (1998) have presented reduced (i.e., con- 
tinuum subtracted) CN images and dust continuum images. These show 
that dust jets and CN jets generally follow the same spiral geometry 
with nucleus rotation. Taken at face value this means that the velocity 
of the CN source and the velocity of the dust is about the same. As 
mentioned already this is in contrast to the case for Halley (A’Hearn et 
al. 1986) where the gas jet and dust jet had different spiral pitch angles, 
indicating different velocities. One interesting exception to the jet coin- 
cidence was a CN-only jet (or possibly more accurately, a gas-only jet) 
reported by Larson and Hergenrother. 

Samarasinha et al. (1998) have presented a preliminary model for the 
CN gas jet which reproduces some of the basic features of their own CN 
images. Their images are otherwise similar to the above-mentioned data 
from Lederer et al. and Larson et al. Although the detailed appearances 
of gas and dust jets change with activity and rotation from day to day 
during this period, the general existence and appearance of dayside dust 
spirals and dayside and nightside gas spirals were consistent. Their 
model has a single active region as the source of the CN jet which 
continues activity on the night side of the nucleus at about 10% of 
the level on the day side. Taken at face value this might be consistent 
with the HCN observations of Wink et al. (1998) mentioned above, 
perhaps suggesting a connection between HCN jets and CN jets, at 
least in comet Hale- Bopp. If the dust jet activity had continued at a 
similar level into the night side, it would be very obvious in continuum 
images. Therefore, it is apparent that the dust jet activity on the night 
side is reduced to well below the 10% level. This is consistent with 
the appearance of dust jets in other comets, even Comet Hyakutake 
(1996B2), which showed well-developed day-side dust jets but fairly 
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uniform (spherical) emission of gas species (CN ar.d OH) from the 
vicinity of the nucleus (Harris et al. 1997). 

Based on the dusty-gas modeling results presented in the previous 
section of this paper, it is clear that the dust-gas dynamic regime 
for comet Hale-Bopp, especially with the dominance of local active 
region(s), is very different from the more typical active comet, e.g., 
P/Halley, and C/Hyakutake (1996B2). Even fairly large particles (10 
/im) in the dust-gas jet will be accelerated to a large traction of the gas 
velocity. This means that it is quite natural for the gas and dust jets 
to be spatially similar in Hale-Bopp, but not so in P /Halley. 

Perhaps the more interesting question raised by the existence of gas 
jets of CN, C 2 , C 3 , NH, OH (implying H 2 0), CO. CS, and H 2 CO, 
and the similar day- night behavior of CN and HCN jets regards the 
original conclusion by A Hearn et al. (1986) that the CN in the CN 
jets of comet P/Halley were produced necessarily from CHON grains. 
As mentioned above, Combi (1987) had suggested that the parent of 
the CN in the jets might have just been a local enhancement of a 
parent molecule, i.e., gas emitted by the local active a*ea which remains 
focused in a jet. A parent molecule would of course move with local 
gas speed, just as the small sub-micron grains invoked to explain the 
CN jets. Dusty-gas hydrodynamic models for localized active regions 
(Kitamura 1986, Korosmezey and Gombosi 1990; Crifo et al. 1995) 
have always indicated that radial expansion near the nucleus is rapid 
and that any jet asymmetry frozen-in by a few nucleus-radii would be 
maintained to considerable distances. Therefore, the general spherical 
symmetry of the gas coma in some comets, is probably not due to a 
simple isotropization from the physics of the inner coma flow expansion 
from day-side active areas. Either distributed produc;ion of gas species 
over the surface or significant extended production of gas from volatile 
grains away from the nucleus is required. Comet Hyakutake (1996B2) 
is an example of a comet with active dayside dust jet production but 
nearly spherical gas (CN and OH) distribution. There is also direct and 
indirect evidence for significant gas production from icy particles away 
from the nucleus (Harris et al. 1997). 

3.3. The Case for C 2 

It has been noticed for a long time that the radial distribution of C 2 rad- 
icals in comets is unusual compared with most other daughter species 
(e.g., CN, C 3 , NH 2 , etc.) whose distributions can all be reproduced by 
a photochemical model that accounts for one generation of exponential 
decay from a parent molecule. Delsemme and Miller (1971) suggest- 
ed that the radial distribution could be explained i: C 2 radicals were 
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trapped in clathrate-hydrate icy grains which were liberated from the 
nucleus similarly to refractory dust particles, dragged away by the out- 
flowing gas, and subsequently vaporized releasing the free C 2 radicals. 
They presented a mathematical icy grain halo model, which could be 
fitted to an observed C 2 profile, as an alternative to the photochemical 
Haser model, which fit poorly. A number of subsequent investigators 
who published spatial profiles of C 2 , sometimes noted the same prob- 
lem and suggested either that C 2 was produced from grains, or it was 
produced as a grand daughter species from a primary parent molecule. 

Despite the difficulty, observed C 2 profiles (Cochran 1985; Com- 
bi &; Delsemme 1986; Wvckoff et al. 1988; Fink, Combi, and DiSanti 
1991; Schulz et al. 1993; Cochran and Trout 1994) and the extraction 
of production rates for all species from standard cometary aperture 
photometry (A Hearn et al. 1995) have typically been analyzed using 
a standard 2-generation Haser model. In recent years more and more 
analyses of observations of OH (e.g., Weaver et al. 1997; Budzien, Fes- 
tou and Feldman 1994; Cochran and Schleicher 1993) and CN (Schulz 
et al.1993 a&b, Schulz et al. 1994) have been performed using a vec- 
torial model. Such in not normally the case for C 2 * The best-fit Haser 
scale lengths have generally depended upon the spatial coverage of the 
data. Parent and daughter scales lengths which are more nearly equal 
are found if the data cover distances close to the nucleus (Wyckoff et 
al. 1988: Fink et al. 1993; Schulz et al. 1993). Parent and daughter scale 
lengths in the ratio of ~6-8 to 1 are found if the data cover distances 
farther from the nucleus (Cochran 1985; Combi and Delsemme 1986). 

By scaling a set of radial profiles of C 2 in comet P/Halley which 
assumed their length scale they follow an r 1 0 variation (r = heliocen- 
tric distance) O’Dell et al. (1988) produced a single average C 2 profile 
covering several orders of magnitude in normalized-distance both very 
close to and very far away from the nucleus. They used a 3-generation 
version of a Haser-type photochemical model with a sequence of three 
scale lengths in the ratio of 1:5:50 corresponding to parent, daughter, 
and C 2 as granddaughter to account for the shape of the composite 
profile. The recent discovery of acetylene (C 2 H 2 ) in comet Hyakutake 
and subsequent measurement in comet Hale-Bopp, provides a prime 
candidate, present at the ~1% level, for the 3-generation chain (C 2 H 2 
to C 2 H to C 2 ). This particular chain has been suggested for many years 
(Jackson 1976; Yamamoto 1981; see the review by Jackson 1991). 

However, this 3-generation model, like the standard Haser model, 
ignores the existence of expected excess energy at each photodissocia- 
tion which has been properly included in the vectorial model (Festou 
1981) and otherwise similar Monte Carlo models (Combi and Delsemme 
1980). In a more recent spatial profile study for individual profiles in 
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comet P/Halley, Combi and Fink (1997) examined ooth 3-generation 
models which included daughter species ejection velocity and a new 
version of a grain source model that modifies some properties of the 
older icy grain halo model (Delsemme and Miller 1971) to be appro- 
priate for CHON grains which have been suggested as sources for CN 
and C 2 jets (A’Hearn et al. 1986, Klavetter and A’Hearn 1994, Schulz 
and A’Hearn 1995). They found that if the 3-generation model is cor- 
rect then there are strict limits to the total ejection velocities for both 
dissociations taken together of < 0.5 km s _1 . In addition, they also 
found that a CHON grain halo could explain the radial distribution for 
a halo size of about 88,000 km at a heliocentric distance of 1 AU. This 
is the expected radius for velocity and lifetime of submicron (~0.2 ftm) 
CHON grains (Lamy and Perrin 1987) if they behave as laboratory 
tholins (Khare et al. 1984). 

Another approach into the question of the source of C 2 radicals has 
been found in comparing the spatial dependence of the (1-0) to the 
(0-0) Swan bands by Roussellot, LaFont and co-workers. Roussellot et 
al. (1994) suggested that their observed variation of C 2 (0-0) and (1-0) 
in comet P/Halley indicates production of vibrationally cool C 2 from 
the nucleus, presumably from the photodissociation of some parent 
molecule and hot C 2 in a spatially extended region and in the vicinity of 
observed jets, which is consistent with production from hot grains. More 
recent observations of comet Hyakutake (C/1996 B2) by Laffont et al. 
(1998) yield similar results indicating the production of cool C 2 near 
the nucleus but the production of high excitation temperature C 2 in the 
region associated with the antisunward-drifting condensations, clumps 
or clusters. Laffont et al. (1998) are presenting the results of similar 
observations and analyses in comet Hale-Bopp at this conference, but 
the results are not yet available as of this writing. 

In an analysis of spatially resolved ID spectroscopic observations of 
Hale- Bopp, Korsun et al. (1998) find similarities between the distribu- 
tions of C 2 and dust continuum not shared by other gaseous species. 
They find that the other gaseous species (NH 2 and CN) are rough- 
ly symmetric (sunward-antisunward) with respect to the nucleus, but 
that C 2 is asymmetric with a sunward preference s milar to the dust 
continuum. They reason that this is due to product on of some of the 
C 2 from the dust. Sekiguchi et al. (1998) also reporl more specifically, 
based on observations of Hyakutake and Hale-Bop ) that most C 2 is 
produced from C 2 H 2 , but that some C 2 is produced from dust. 
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3.4. Dust Fading and Fragmentation 

A radial profile of column abundance (or brightness in the optically thin 
limit) varying as the inverse of the projected distance, p, is found for a 
simple constant point source expanding at a constant velocity. This is 
the simplest picture of dust production in comets. Even in the realistic 
case where solar radiation pressure is important, the radial profile in the 
direction away from the sun, as projected on the sky plane, is still found 
to vary as the inverse of the projected distance (Wallace and Miller 
1958). Evidence for the departure from a 1/p distribution were reported 
by Jewitt and Meech (1987) and Jewitt and Luu (1989) from analyses 
of the radial profiles of dust continuum observations in comets. In spite 
of a technical error in their mathematical model discussed by Combi 
(1994), the finding by Baum et al. (1992) of what they called dust fading 
on scales of ~ 10 4 - 10 5 km in some comets, but not all, appears to 
remain valid. They found departures from the simple 1/p distribution 
which could be interpreted as a loss of particles with increasing distance 
from the nucleus. Preliminary indications from observations of comet 
Hale-Bopp (Schleicher, private communication) are that grain fading 
does not occur on this spatial scale. 

On the other hand, there is evidence that grain fragmentation in 
the very innermost comae of active comets may be quite important. 
The brightness profiles of dust continuum determined with the Giotto 
Halley Multicolor Camera experiment (Keller et al. 1986) within a few 
tens of km of the nucleus of comet P/Halley flattened considerably. This 
occurred at distances where dust particles could still be accelerated by 
outflowing gas, and one might expect to have observed the brightness 
profile become more steep not less. The fragmentation of particles has 
the opposite effect from fading on the radial profile (Thomas and Keller 
1990; Konno et al. 1993). As particles break apart, the total surface 
area of all the small particles is larger than those of the original large 
particles, in the absence of evaporative loss. This increases the total 
scattering cross section and thus increases the brightness, relative to 
the base-line 1/p with increasing distance. 

In addition to the Halley Giotto images, there is evidence for frag- 
mentation in the analysis of ground-based dust coma and tail images. 
In an analysis of dust tail images of comet Austin (1990V) with their 
inverse Monte Carlo model, Fulle et al. (1993) found that a considerable 
spread of initial velocities, the so- called terminal velocity, was required 
for particles having any particular size. This is contrary to the standard 
Finson-Probstein (1968) dusty-gas dynamics picture which implicitly 
assumes non-fragmenting particles, and where particles of a certain size 
and bulk density are accelerated by the gas flow to a particular single 
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velocity, i.e., the velocity distribution for a certain size particle is a delta 
function. Combi (1994), in order to reproduce ground-based dust coma 
images of comet Halley using a Finson-Probstein-t/pe Monte Carlo 
simulation, introduced a heuristic dust-fragmentatkm function, which 
assumed that particles started with some velocity given by dusty-gas 
dynamics and then fragmented into a population of smaller particles. 
More recent modeling by Fulle et al. (1997) for Hyakutake (1996B2) 
similarly requires a velocity dispersion for each particle size, implying 
that particle fragmentation is important within, and perhaps just out- 
side, the dusty- gas acceleration region (typically a couple of hundred 
km). 

Based on the modeling discussed in the previous section of this 
paper, it appears that for comet Hale-Bopp a wide range of particle 
sizes will be accelerated to nearly the gas velocity, especially within the 
concentrated jets. The more uniform velocities for the particles in the 
0.1 to 10 /i m range, which typically dominate the observed visible dust 
coma, will probably make the analysis of dust coma and tail images a 
less sensitive indicator of particle fragmentation than other recent and 
more moderate comets. All particles regardless of whether they were 
primary parent particles, or fragments of them, will all be accelerated 
to nearly the gas velocity and thus will not have much dispersion in 
velocity. 


4. Summary 

We have presented a review of the state of knowledge legarding issues of 
dust-gas interrelations in comets, especially as it relates to comet Hale- 
Bopp at the time of this writing. We include (1) theoretical modeling 
of the gas/dust environment for Hale-Bopp, (2) a review of certain 
important relevant issues as posed from the results of comets observed 
before Hale-Bopp, and (3) a discussion of new Hale-Bopp observations 
which either shed new light or at least herald the passible shedding of 
new light in the not-too-distant future. Hale-Bopp is still being actively 
observed by the world l s major observing programs. Despite the fact 
that we now only have an incomplete and preliminary picture of the 
observational evidence in this important area, there are already many 
interesting and important findings and questions. 

The dust-gas environment for comet Hale-Bopp especially when 
very active (r < 1.5 AU), is very different from the typical bright comets 
(e.g., P/Halley, and C/Hyakutake 1996B2). The large dust-gas produc- 
tion has a number of important consequences. 

- The collisional fluid region is large. 
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- The resulting increased photochemical heating efficiency signifi- 
cantly increases gas outflow speeds and their variations with distance 
from the nucleus and with heliocentric distance. 

- The large collision region will have an effect on the applicability 
of free-molecular- flow vectorial models which implicitly assume colli- 
sionless expansion. 

- The UV opacity in the inner coma is not negligible, and could 
effect parent lifetimes and ion production especially for observations 
made with high spatial resolution. 

- The high UV opacity will shield ionization within about 3000 km 
from the nucleus and within a tailward shadow (at r = 1 AU), and alter 
typical ion-neutral chemistry. 

- The dust-gas interaction is strong and the region where it is impor- 
tant is large. 

- This interaction yields fast micron-sized dust particles accelerated 
to nearly the gas speed. Even larger mm-size particles are accelerated 
to much larger velocities than the typical bright comet. Such effects 
will be important in understanding dust kinematics and in determin- 
ing total dust-to-gas mass ratios from long wavelength IR and radio 
observations. 

Preliminary observational results regarding dust-gas interrelations 
are also very interesting. 

- The production of gas species and dust seems to be strongly influ- 
enced and sometimes dominated by localized active regions that rotate 
with the nucleus. 

- The dust production from the active regions essentially turns-off 
on the night side. 

- The gas production is reduced on the night side, but not completely 
inactive (~10% level in some species, but CO is uniformly active). 

- There is evidence for extended sources not only for CO and H 2 CO, 
as seen before, but also for CH 3 OH, and HNC. 

- At large distances (r > 4 AU) cometary CN seems to be mostly 
produced from HCN. At smaller distances (r < 4 AU) another source 
(or sources) appears to become active; suggestions for grains or other 
parent molecules (C 2 N 2 , CH 3 CN) have been made. 

- There is evidence for multiple sources for cometary C 2 : photodisso- 
ciation of C 2 H 2 is an obvious candidate, and production of vibrationally 
excited C 2 from grains has been suggested. 

- The new list of gas-jet species from visual-range CCD images of 
comet Hale-Bopp includes not only CN, C 3 , and C 3 as in Halley, but 
also NH, and OH. Doppler-resolved radio observations indicate a simi- 
lar behavior for CO, H 2 CO, CH 3 OH and H 2 S. Similarities of behavior 
of HCN jets in Doppler-resolved radio observations and CN jets in CCD 
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images, and the presence of gas jets in most commonly observed species, 
especially OH and NH, open to question the conclusion that previous- 
ly observed CN jets necessarily implied production from an extended 
source of CHON grains. We examine here the pros and cons of three 
scenarios: 

1. Gas-jet species are produced from CHON grains. Pros: Small 
organic grains (sub-/xm) would be accelerated to nearly the gas speeds 
in both Halley and Hale-Bopp, accounting for the similar spatial geom- 
etry of dust and gas jets in Hale-Bopp but also their dissimilarity in 
Halley. Small organic grains (0.2 mm) should live long enough (1-2 
xlO 5 seconds) for a comet at 1 AU (Lamy and Perrin 1988) to provide 
a spatially extended source. Cons: The new evidence of gas jets in OH 
and NH in Hale-Bopp indicates the gas-jet source may be rooted in 
the volatile fraction of the comet rather than in organic CHON grains. 
The similar visibility of OH jets and CN jets indicates that the relative 
source rates are similar in the jets source as in the comet in general. 
This similar relative abundance was also the case for the OH, CN and 
C 2 arc-structure seen in the tail of comet Hyakutake (C/1996 B2) which 
indicated production by typical fragments of the bulk nucleus rather 
than just CHON grains which are only one component of cometary 
material. 

2. Gas-jet species are produced from volatile nucl eus fragments (icy 
grains). Pros: This picture accounts for production of OH and NH 
which are considered to be produced from the volatile fraction of the 
cometary nucleus material, rather than from the organic CHON grains. 
A similar process has been invoked to explain the production of OH, 
in addition to CN and C 2 from a secondary source cown in the tail of 
comet Hyakutake (C/1996B2). Cons: In Hyakutake, both the lifetime 
of the particles (from the persistence of the source) ind the motion of 
the source down the tail indicated that the parent ‘ragment particles 
were large (mm-sized). In Hale-Bopp the gas-jet scurce particles are 
probably still closer to fim size rather than mm size based on the fast 
velocity and the fact that small icy grains do not survive more than 
about 10000 km from the nucleus (Delsemme and Miller 1971). Of 
course, the particles responsible for the jets in Hale-Bopp do not have 
to be the same as those responsible for the tail-side arc in Hyakutake. 

3. Gas-jets axe produced simply from localized con :entrations of par- 
ent molecules. Pros: Because parent molecules are p art of the general 
gas outflow, they, like the smallest grains, will move with the general 
gas outflow velocity. This would account for the velocity sorting for 
jets in Halley and for the lack of it in Hale-Bopp. The dusty-gas hydro- 
dynamic models indicate that non- spherical asymmetries in the coma 
are frozen-in by a few nucleus-radii from the nucleus, therefore con- 
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centrations in angle of parent molecules could remain concentrated to 
large distances, thereby possibly maintaining a gas-jet identity. Cons: 
The gas concentrations which are similar to dust concentrations near 
the nucleus (10° half-cone-angle FWHM Gaussians) tend to broaden 
to ~30° half- cone angle Gaussians by a few nucleus radii. Are these 
narrow enough to account for the well focused gas jets? 

- A few more general and difficult questions regarding both gas and 
dust jets remain. It may be difficult to account for the constant project- 
ed radial propagation velocity of gas and dust jets (0.4 km s -1 ) solely 
on geometrical projection effects (Larson, private communication). This 
has larger consequences than the exact nature of the gas and dust jets. 
It is clear from the modeling presented here that the outflow velocities 
(0.7 to 2 km s" 1 ) predicted by dusty-gas dynamics models are con- 
sistent with those observed in Doppler resolved measurements (Biver 
et al. 1998; Colom et al. 1998). How can gas jets propagate at much 
smaller speeds than the actual gas outflow speed? 

Over the coming couple of years the wealth of important observa- 
tional results for comet Hale-Bopp will be filling the scientific litera- 
ture. It is very likely that as a result of this singularly bright comet 
we will obtain definite answers to many of the above important ques- 
tions about dust-gas interrelations. They have important consequences, 
not just those pertaining directly to our understanding of the behavior 
of processes occurring in the coma, but also, and more significantly, 
on our understanding of the detailed compositional and morphological 
make-up of comets, and how and where they were formed. 
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Figure 3. folate 1. 3D Dusty-Gasdynamic Model of a Dust-Gas Jet. Shown are (a) 
the nilhtber density of gas and (b) the speed of the gas. The calculation corresponds 
to comet. Hale-Bopp at ^1 ALL The jet is directed toward the right. 
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Figure Plate 2. 3D Dusty-Gasiivnamic Model of a Dust-Gas Jet. Shown are (a) 
the bulk mass density and (h) the speed for the 1 |un dust particles. The calculation 
corresponds to comet Hale-B«>pp at AU. The jet is directed toward the right. 
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Figure 5. Plate 3. 3D Dusty - G asdy namic Model of a Dust-Gas Jet. Shown are (a) 
the bulk mass density and (b) the speed for the 1 mm dust particles. The calculation 
corresponds to comet Hale-Bopp at -^1 Al T . The jet is directed toward the right. 


DusGasl0.tex; 15/09/1998; 14:37; no v. ; p.33 
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Figui'e 6. Plate 4. 3D Dusty -Oasdynamic Model of a Strong Dust-Gas Jet. Shown 
are (a) the number density of t he gas and (b) the bulk mass density for optical t 
/i m particles in the vicinity of a strong jet, directed toward the right. 
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Abstract. The first results for applying a three-dimensional multiscale ideal MHD model for the 
mass-loaded flow of Jupiter’s corotating magnetospheric plasma past Io are presented. The model 
is able to consider simultaneously physically realistic conditions for ion mass loading, ion-neutral 
drag, and intrinsic magnetic field in a full global calculation without imposing artificial 
dissipation. Io is modeled with an extended neutral atmosphere which loads the corotating 
plasma torus flow with mass, momentum, and energy. The governing equations are solved using 
adaptive mesh refinement on an unstructured Cartesian grid using an upwind scheme for MHD. 

For the work described in this paper we explored a range of models without an intrinsic magnetic 
field for Io. We compare our results with particle and field measurements made during the 
December 7, 1995, flyby of Io, as published by the Galileo Orbiter experiment teams. For two 
extreme cases of lower boundary conditions at Io, our model can quantitatively explain the 
variation of density along the spacecraft trajectory and can reproduce the general appearance of the 
variations of magnetic field and ion pressure and temperature. The net fresh ion mass-loading rates 
are in the range of —300-650 kg s~ \ and equivalent charge exchange mass-loading rates are in the 
range -540-1 150 kg s' 1 in the vicinity of Io. 


1. Introduction 

Io s interaction with Jupiter’s corotating inner 
magnetosphere has been studied for over 25 years. Io has a 
neutral atmosphere which is probably locally thick but rather 
uneven across its surface [Lellouch et al, 1992; Ballester et 
al - , 1994]. The ultimate source for atmospheric gases appears 
to be the numerous active volcanoes on the surface. It is the 
energetic particle environment near Io which is responsible for 
the balance of the plasma heating. Joule heating, ionization, 
and surface and atmospheric sputtering, and which in some 
form drives the escape of the neutral atmosphere [Schneider et 
qI., 1989]. A realistic model for the flow of magnetospheric 
plasma, interacting with lo’s neutral atmosphere is necessary 
for a complete understanding of the global picture. 

The results of the measurements by the particle and field 
instruments on the Galileo Orbiter during the December 1995 
flyby of Io provide new and important information with w hich 
a realistic MHD simulation for the plasma interaction can be 
tested. Kivelson et al [I996a,b] have recently suggested that 
an intrinsic magnetic field for Io with a strength and direction 
which would nearly cancel the local Jovian field near the 
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surface of Io, can explain many aspects of the Galileo 
magnetometer measurements taken during the December 1995 
flyby of Io. However, the large plasma densities and low 
temperatures measured at the same time by the Galileo plasma 
science package [Frank et al., 1996] and the plasma wave 
instrument [Gurnett et al, 1996], indicate that ion mass 
loading is significant and may provide an explanation for the 
magnetic field measurements without requiring an intrinsic 
field. 

Our understanding of the interaction of Io with Jupiter’s 
magnetosphere goes back to the discovery of Io-related 
decametric radio emission discovered by Bigg [1964] and the 
unipole inductor model which explains it [Goldreich and 
Lynden-Bell 1969]. There have been a number of theoretical 
studies of this interaction during the immediate post-Voyager 
era [e.g., Southwood et al, 1980; Neubauer, 1980] (see also 
Hill et al [1983] for an excellent post-Voyager review). 
Neubauer [1980] presented an analytical model of Alfven 
standing wave current system which connects current through 
the ionosphere of Io. Southwood et al [1980] examined data 
from several Voyager instruments and examined the possible 
role of an intrinsic magnetic field for Io as a way to retain a 
robust enough ionosphere to provide enough conductivity for 
completing the Io-Jupiter current circuit. Early theoretical 
work was often done either in the context of a “thin” 
atmosphere [e.g., see Cloutier et al, 1978] indicative of the 
surface temperature (-130 K), or “thick” extended neutral 
atmosphere [e.g., see Goertz, 1980] more indicative of volcanic 
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temperatures (-1000 K). Subsequent evidence [Sieveka and 
Johnson , 1985; Schneider et ai, 1991; Leltouch et ai, 1992; 
Ballester et ai, 1994; Belton et al., 1996] seems to indicate a 
mixed picture of the global atmosphere, which has a large 
extended corona like a thick atmosphere, but appears to be 
dominated by local major injection of hot (high speed) 
gas/dust plumes to high altitudes but only near active 
volcanic vents. Therefore, although the atmosphere is 
probably only locally thick, it still has a large extended 
neutral corona which might provide a sufficient source of 
impact ionization and photoionization. 

There have been a number of three-dimensional (3-D) 
numerical studies of the plasma flow past Io [Linker et ai, 
1988, 1989, 1991; WolfGladrow et al., 1987]. Because of the 
complexities in the problem and numerical difficulties, studies 
have been limited to restricted portions of the problem, have 
introduced fixes such as artificial viscosity, or have been run 
using less than fully physically realistic conditions. Linker et 
ai [1991] have performed the best example to date of 3-D 
MHD modeling of lo. In this model they used the two-step 
Lax-WendrofT method. These calculations were performed on a 
nonuniform but fixed structured grid, and artificial dissipative 
terms (artificial viscosity) had to be introduced to the basic 
MHD equations in order to make them computational stable. 

A new numerical method has been developed at the 
University of Michigan for solving the full 3-D MHD 
equations which resolves a number of difficulties that have 
plagued past attempts. The 3-D version of this new method, 
called multiscale adaptive upwind scheme for MHD (MAUS- 
MHD), has already been applied to the full comet-solar wind 
interaction problem and has been discussed in detail by 
Gombosi et ai [1996]. This interaction is quite similar to that 
expected for the Io-Jupiter-magnetosphere interaction in the 
absence of an intrinsic magnetic field for Io. A full Earth 
magnetosphere version of the same calculation has also 
recently been developed [De Zeeuw et ai, 1996], so the method 
is also able to treat intrinsic magnetic fields as well, if 
necessary. We can and have run test cases with an intrinsic 
field for Io, but in this paper, concentrate on trying to explain 
the existing data with nonintrinsic field models. 


2. Model Description 


hydrodynami :s; MAUS-MHD, manuscript in preparation, 
1998]. In th s scheme div B/r is treated as a passive scalar so 
that any numerically generated nonzero contributions are 
convected a»vay. Once div B = 0 is imposed as an initial 
condition tc the problem it is satisfied to within truncation 
error. 

Our model of Io’s interaction with Jupiter’s magnetosphere 
is essentially similar to the comet problem of Gombosi et ai 
[1996]. The dimensionless conservative form of the ideal MHD 
equations can be written as 



where W an d S are eight-dimensional state and source vectors, 
while F is an 8 x 3 dimensional flux diad. All quantities 
denoted by i tilde are normalized by physical quantities in the 
undisturbed upstream region 


t — a<x>t / R Jq 

(2a) 

? = r/R Io 

(2b) 

p = p/p- 

(2c) 

u = u/a„ 

(2d) 

p = p/p»a„ 2 

(2e) 

B — B / ^jiQpooaoo 2 

(20 

where t is lime, r is radius vector, p is mass density, 

u is bulk 


flow velocity, p is pressure, is the permeability of free space, 
and B is magnetic field vector. The subscript « refers to values 
in the undisturbed upstream flow. Finally, Rj 0 is the radius of 
Io and aoo is the hydrodynamic sound speed in the undisturbed 
upstream pi isma flow. 

The normalized state and flux vectors are written as 


In the 3-D MAUS-MHD model the calculation is done on 
an unstructured Cartesian grid using octree adaptive mesh 
refinement (AMR). This makes optimal use of computer 

resources allowing problems with disparate scales and 

unusual and varying geometries to be modeled globally, while 
still retaining sufficient resolution locally to capture shocks 
and other important features resulting from naturally occurring 
sharp gradients in the flow. The octree structure is a 
hierarchical cell structure, based on multiple generations of 
cubic parent cells which can be split into eight daughter cells. 
The adaptive data structure has been explained by DeZeeuw 
and Powell [1992]. The approximate eight- wave Riemann 
solver for ideal MHD is combined with a second-order 
MUSCL-type scheme [van Leer , 1979]. It uses a novel 
approach to handle the “div B problem” derived and 
discussed at length by Powell [1994] and Powell et ai [1995, 
also Multiscale adaptive upwind scheme for magneto- 
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where £ is the normalized internal energy density 

pu 


- 1 
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(4) 


The soiree vector in (1) corresponding to ion mass-loading, 
M, is of ti e same form as used by Gombosi et ai [1996], and is 
given in normalized units as 
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Qi is the ionization rate per unit volume evaluated at Io’s 
radius; F in is the ion-neutral collision mo men turn -energy 
transfer rate per unit volume, also evaluated at Io’s radius. In 
the simplest interpretation Q,(R [o /r) a is given by the neutral 
density, n n , divided by x\ y the ionization timescale for neutrals, 
F in (Rio/r) a is given by n n kj n , where kj n is the ion-neutral 
momentum transfer collision rate, and n n is the neutral density. 
The choice for the exponent, a, is discussed below. For the Io 
cases presented here the net neutral flow velocity, u n , is taken 
to be zero. This could be generalized in the future to account 
for a net escaping neutral atmosphere at high altitudes. 

For the Io version of the MAUS-MHD model simulations 
were performed both for a pure conducting sphere with and 
without ion mass loading and for a sphere with fixed boundary 
conditions, as discussed in detail below. The simulation 
volume is 900 x 600 x 600 R Io , and there are 9 levels of 
refinement and 92,000 cells with sizes ranging from <0.1 R !o 
near Io to 50 R Io far upstream and downstream. We have 
chosen the upstream plasma conditions as follows from our 
best estimates of the Galileo December 1995 epoch, including 
the various recently published Galileo measurements (see 
Table 1). Since the calculation is performed in normalized 
units, there is some flexibility to rescale the results over a 
range of realistic conditions, but this cascades to other 
parameters [see, e.g., Gombosi et ai, 1996]. We have 

investigated numerous versions of cases for a simple 
conducting sphere, a mass-loading ionosphere, and with an 
intrinsic magnetic field for Io. 

The formulation of mass loading for comets [ Gombosi et al., 
1996] permits us to incorporate both the effects of mass loading 


Table 1. Upstream Plasma Conditions 


Parameter 

Value 

Upstream plasma density 

3500 cm' 3 

Upstream plasma temperature 

92 eV 

Upstream mean molecular mass 

22 amu 

Upstream magnetic field 

1800 nT 

Corotation flow speed 

56.8 km s ! 

AlfV^nic Mach number 

0.4 

Mach number 

2.2 

Ratio of specific heats (y) 

1.667 


from new ionization (impact and photoionization) as well as 
momentum and energy loading from charge exchange and 
nonreactive ion-neutral collisions. In the comet calculation 
Gombosi et al. used the term “ion-neutral friction” to describe 
the effect of collisions which cause momentum and energy 
exchange but which result in no new ions. The dominant 
source of such an exchange of momentum and energy results 
from symmetric charge exchange reactions involving the 
dominant torus ions and atmospheric neutrals (O^ + O, S + + S, 
and + S). Linker et al. [1991] appropriately included an 
ion mass loading rate and a “charge exchange rate" in their 
MHD models. However, many nonreactive collisions between 
ions and atmospheric neutrals also can exchange momentum 
and energy without yielding net new ion mass. These so- 
called knock-on collisions are important for producing the 
extended neutral corona of Io [Sieveka and Johnson , 1985; 
Johnson, 1990]. Furthermore, it has been shown by Smyth and 
Combi [1997] that the distribution of sodium in the corona 
[Schneider et ai, 1991] as well as the more familiar B-cloud is 
a natural outcome of the collisional cascade process 
(incomplete atmospheric sputtering) that is initiated by such 
collisions. For simplicity in the remainder of this paper this 
momentum-energy loading term is referred to as the charge 
exchange rate. 

On the basis of Voyager-era conditions, Smyth [1992] has 
estimated a total collision rate of -2 x 10 2 * s' 1 between 
incident plasma torus ions and neutrals in lo's atmosphere 
with about 1/3 yielding fresh ions deposited into the plasma 
flow and the rest contributing to momentum/energy exchange 
alone, i.e., the charge exchange rate. These estimates are based 
on extrapolations of knowledge gained about the extended 
sodium atmosphere of Io to the more abundant neutral species 
S, O, SO, and S0 2 . On the basis of their own estimates. Linker 
et al. [1989, 1991] have explored MHD models with a 

comparable charge-exchange rate but with a fresh ion mass- 
loading rate nearly an order of magnitude smaller (~10 2? s' 1 ) 
than Smyth’s. Based on the then non-detection of optical/tJV 
emissions close to Io, Shemansky [1980] provided an upper 
limit to the fresh ion mass-loading rate of - IO 27 s* 1 . Smuh and 
Shemansky [1983] later raised this to 4 x IO 27 s‘ l after 
subsequent detections of neutral O and S emissions from 
ground-based and space-based remote measurements. It is also 
likely that the torus is more robust in recent years, and that 
mass, momentum, and energy input rates may all be higher now 
than 17 years ago. This appears to be borne out bv 
comparisons between Galileo measurements during the 
December 1995 flyby and Voyager measurements [Frank et al , 
1996; Gurnett et al . , 1996]. 

As is evidenced by the range of values which appear in the 
literature, making an independent estimate for ion mass loading 
and charge exchange rates is difficult, because they are based 
on a number of uncertain and indirect measurements and 
assumptions. With the availability of the Galileo PI S 
measurements [Frank et al., 1996] of plasma densitv 
temperature, and pressure we have adopted the different 
approach of using the net ion mass-loading and charge 
exchange rates as adjustable parameters and constraining the 
values based on comparison of our model results with the 
measurements along the spacecraft track. 

S0 2 appears to be the major primary constituent in lo's 
atmosphere [Ballester et al, 1994; Lellouch et al., 1992], but 
observations still to date only provide hemispherical l> 
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averaged total column abundances. The sodium eclipse 
measurements of Schneider et ai [1991] are the only 
unambiguous information we have about the vertical 
distribution of neutral gas species in the corona of Io from 
about 1.4 to about 6 Rj 0 . The interpretation of recent Hubble 
observations of various UV emissions of atomic S and O 
[Bailester et ai, 1996] is complicated by the fact that the 
spatial distributions of the emissions are determined by a 
combination of the mass distribution in the atmosphere with 
excitation by the plasma density and temperature environment. 
Theoretical and empirical models have explored the structure 
of the neutral atmosphere and its behavior to various gas 
sources and energetic heating and ejection processes, such as 
sublimation, volcanic ejection, surface and atmospheric 
sputtering, UV, plasma, and Joule heating, and LTE and non- 
LTE IR radiative cooling [Sieveka and Johnson , 1985; 

McGrath and Johnson , 1987; Moreno et ai . 1991; Strobel et 
ai, 1994; Pospiezalska and Johnson, 1996; Wong and 
Johnson , 1995; Smyth and Combi, 1997 ]. Despite recent 
advances in observational and theoretical work, still only a 
rough sketch of the structure of the neutral atmosphere can be 
made. 

For our MHD simulations we have adopted ion mass, 
momentum, and energy loading source terms in the vicinity of 
Io as if proportional to the neutral density distribution which 
is scaled as a power law in distance (r/R [o ) with an exponent, a. 
We do not need to specify the neutral density itself, however, 
just the ionization mass and momentum/energy exchange rates 
and their dependence on distance. We have adopted the power 
law with an exponent of a - -3.5 as found by Schneider et ai 
[1991] for the trace species, sodium, outside about 1.4 R Io . 
Smyth and Combi [1997] have shown that this vertical 
distribution of sodium from the Schneider et ai. [1991] eclipse 
observations and the more typical emission profiles of more 
spatially extended sodium can together be explained by a 
single modified incomplete cascade sputtering distribution, 
which results naturally from multiple neutral -neutral and ion- 
neutral collisions between torus ions and atmospheric 

neutrals. Therefore it is reasonable to expect a similar type 
distribution for neutral species in general. This power law is 
contrasted with the less steep inverse square-type power law 
distribution assumed in past MHD mass-loading calculations 
[Linker et ai, 1988, 1991]. For the mass-loading models 
presented here, we assume the neutral distribution has this 
dependence from an inner boundary at 150 km above the 
surface of Io, out to the Lagrange radius of 6 R [o . 

If details about the ion composition and electron 
temperature in the torus, and the global spatial distribution 
and composition of the neutral atmosphere were available, it 
would be possible to make a better estimate of the spatial 
distribution of the ion mass-loading rate and the charge- 
exchange rate (the collisional energy-momentum loading rate) 
needed to specify the source terms in the MHD equations. 
Typical cross sections [McGrath and Johnson , 1989] between 
oxygen and sulfur ions dominating the torus and O, S, and SO 2 
neutrals dominating the atmosphere 'how that typical 
symmetric charge exchange reactions have cross sections in the 
range of 20 to 40 A 2 whereas most impact ionization and 
elastic collisions [Johnson, 1990] have cross sections in the 
range of 1 to 20 A 2 . Electron impact of course only nets ion 
mass but is highly dependent on the electron temperature. 

The power law distribution may not be correct well inside 
of 1.4 R Io because the range of the observ ations of Schneider et 


ai [1991] started at this distance. At much lower altitudes the 
resulting ion mass-loading rate is a trade off between the likely 
rapid increasing of neutral density (probably more like an 
exponentia scale height distribution), with the decreasing 
plasma flu < (electron and ions) caused by recombination and 
plasma-net tral collisions themselves. Flctser et ai [1997] 
have shown preliminary results of Galileo radio occultation 
measurements which indicate electron densities in excess of 
50,000 cm 3 at the peak of the ionosphere, -50-80 km above the 
surface. This is similar to but higher than given by Kiiore et 
ai [1974, 1975] from the Pioneer era, and again generally 
consistent with the more robust torus of the recent times. At 
this time, a fully coupled MHD- ionosphere calculation would 
be beset with many assumptions, for example, the neutral 
atmosphere density profile which is still uncertain to within a 
factor of a few. In any case the scope of this type of study is 
currently oeyond the state of the art even for the Earth for 
which there is much more knowledge. 

For the models presented in this paper we have adopted two 
types of lower boundary conditions which represent opposite 
pictures of the type of conditions present above Io’s 
ionosphere One is a perfectly conducting impenetrable 
sphere. In this case both the magnetic field and the plasma flow 
are reflected from the lower boundary. This represents an 
extreme case of not allowing either plasma flow or magnetic flux 
to penetrate the sphere, and we hereafter refer to this as a 
reflecting boundary condition. The second case corresponds 
to fixed boundary conditions where the magnetic field and 
plasma density and velocity are specified in ghost cells just 
below the boundary. In this case both magnetic flux and 
plasma flow can pass through the boundary. For the fixed 
boundary condition model presented here the plasma density 
is set to a value of 3 times the upstream torus density (or 
10,500 cm' 3 ) which is representative of the global spherical 
average indicated by radio occultation measurements at an 
altitude ol 150 km. The plasma velocity is set to zero, being 
consistent with strong coupling to the neutral atmosphere. 
The choice of a value for the fixed magnetic field at the lower 
boundary requires some discussion. We hereafter refer to this 
case as fixed boundary conditions. 

Linker u ai [1991] presented model results which solved 
for the full MHD equations (mass, momentum, energy’, and B 
field) outsi de of Io, and only the magnetic field with a constant 
conductivity inside the inner boundary, throughout the 
interior of lo. The internal conductivity was set in such a way 
to obtain a high enough integrated conductivity (10-150 
mhos) neci ssary for maintaining the Jupiter-Io current system. 
They point out that it is not really known where the current 
loops clos j, i.e., in the ionosphere but outside of Io, or even 
through the dense lower ionosphere and into a possible 
conducting interior. Linker et al. found, however, that the 
same exte ior solutions could be obtained whether they 
performed the computationally intensive calculation of the 
magnetic f eld inside Io or whether they simply held the B field 
at the inne ■ boundary constant at the background level. In fact, 
their full calculation showed that in any case “the magnetic 
field changed very little inside the sphere.” 

The physical interpretation for a constant magnetic field is 
that if Io does not have an intrinsic B field (i.e., from an internal 
dynamo), t would be expected that a value similar to the 
average in posed Jovian external field would be induced over 
time. We lave estimated that the time-averaged value of the 
field at Ic happens to be slightly less (-98%) than the 
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instantaneous value at lo's position at the point of Galileo s 
close approach. The same ratio was found either roughly using 
the variation from the old Goertz el al. [1976] magnetic field 
model or more carefully using the newest model of Connerney 
el al [New model of Jupiter’s magnetic field constrained by 
the lo flux tube footprint, manuscript in preparation, 1998]. It 
is noteworthy that those models predict values that are overall 
too low and too high, respectively, compared with the actual 
values measured by Galileo [Kivelson et at.. 1996b). 

As mentioned already for this study we have left the ion 
mass-loading rate and the charge exchange rate (pure 
momentum/energy exchange) as adjustable parameters. We 
then undertook an iterative procedure to vary these quantities 
and compare our model results with the measured plasma 
density, pressure, and magnetic field along the Galileo 
spacecraft trajectory until a reasonable (albeit not perfect) 
match was obtained. Upstream we assume the aforementioned 
torus conditions. Downstream and to the sides of the 
simulation box we assume absorbing boundary conditions, i.e., 
mass, momentum, and energy, exit, and nothing comes back. 
Because of the adaptive grid technology we are able to place 
the upstream, downstream, and lateral boundaries far enough 
from lo so they do not affect on the solution. 


3. Model Results 


To serve as a baseline, we show results for a perfectly 
conducting sphere without mass loading, given realistic 
upstream plasma torus conditions. We refer to this case as no 
mass loading and reflective boundary conditions. This is 
shown mainly for comparison with the similar nonintrinsic 
field model of Kivelson et al. [1996b], and to serve as a point of 
comparison with the mass-loaded cases. For this case, the only 
parameters to be set are the upstream Mach number of 2.2 and 
the upstream Alvdnic Mach number of 0.4. In (5) the ion-mass 
and ion-neutral drag rates Q, and F,„ are zero. The results are 
scaleable to other values of upstream magnetic field and 
conducting radius given constraints imposed by the Mach 
number (which does not have a great effect on the overall 
global structure of the flow) and the Alfvemc Mach number 
(which is more critical). In any case our estimates of both the 
upstream Mach number and Alfvenic Mach number are 
consistent with the measured Galileo flyby conditions and 
rigid corotation. 

We show results for two mass-loading cases where the 
particular values of the ion mass-loading rate and the charge- 
exchange rate serve as the adjustable parameters. We have run 
some models with larger exponents (-4 to -5), however, -3.5 is 
based on an actual measurements and is quite satisfactory. he 
goal is to be able to reasonably reproduce the variation of the 
plasma density, pressure (and temperature), and magnetic field 
along the spacecraft trajectory according to the results of the 
Galileo lo flyby results [Frank el al . 1996; Kivelson et al. 
1996b; Gurnet I et al. 1996], The sensitive features in the 
plasma data to be reproduced are ( 1 ) the height of the plasma 
density peak along the center of the wake, (2) the reduction of 
the magnetic field, (3) total width of the disturbance m all 
quantities, (4) the reversal and side peaks ot the pressure and 
temperature, and (5) the vector velocity magnitude and 
direction. Frank et al. [1996] have reported plasma velocities 
(44 km s' 1 ) that are below corotational speeds (56.8 km s ) 
even 4 Rjo from the central axis of the wake. However, these 


have a large enough uncertainty (-5 and +15 km s ’) to be 
consistent with corotational speed, which seems more 
reasonable. 

Plates la-lf show a comparison of slices of the magnetic 
field and flow streamlines in the planes parallel and 
perpendicular to the upstream magnetic field for the non-mass- 
loaded and both mass-loaded cases. In the non-mass-loaded 
reflective case one can see the strong influence of the Alfv£n 
wings in the parallel plane (la) [see Neubauer. 1980; Linker et 
al 1991] as the flow diverts around the conducting sphere 
and as the disturbance is directed off at the expected angle (tan 
9 = 1/M A ) where 0 is the deflection angle and M A is the 
Alfv6nic Mach number. In MHD the Alfv6n wings are the 
analog of vortex sheets of regular fluid dynamics [Kogan, 
1961]. In the parallel plane for the mass-loading case with 
reflective boundary conditions (lb), the effect of the inner 
boundary conducting sphere, which dominates the non-mass- 
loaded case, is weakened but still present. Finally, in the 
fixed-boundary mass-loading case (lc) the effect of the Alfven 
wings has all but disappeared. Plates Id- If show the 
comparable plots in the plane perpendicular to the upstream 
magnetic field. The December 1995 Galileo flyby trajectory lay 
very close to this plane. The diverted flow streamlines are not 
very different between any of the three cases, in fact for the two 
mass-loaded cases the flow fields in this plane are virtually 
identical The small variations in the magnetic field magnitude 
are apparent. Note that in mass-loading fixed case (Plates lc 
and If) the diversion of flow is similar in both planes. 

Plates 2a, 2b & 2c show the striking differences for the 
plasma ion density among the various cases. The dark lines 
show the outlines of the octree cell structure of the adaptively 
refined mesh. In the non-mass-loading reflective case (2a) the 
flow diversion around the sphere is evident. The density 
concentrates somewhat along the Alfvdn wings but the wake is 
essentially evacuated of material. The mass-loaded reflective 
case (2b) shows a large density pileup in the wake right 
behind lo and a fair dispersion of material spread in the 
direction parallel to the upstream magnetic field. This is 
caused by the strong influence of the reflective boundary 
conditions. The mass-loaded fixed boundary case (2c) also 
produces a plasma wake, but it is more concentrated along the 
central line of the wake, does not have the large extension 
along the field-parallel direction (as in Plate 2b), and is more 
extended along the wake. 

Direct comparisons of the model with the Galileo t > y 
data, along the trajectory are given in Figures la- Id. Model 
results are superimposed on the measurements for both mass- 
loading models with reflective and fixed boundary conditions. 
The reduction of the magnetic field in the wake along the 
spacecraft trajectory is shown in Figure la. Both of these 
models provide a similar breath and depth of the magnetic field 
disturbance as the pure intrinsic dipole model of Kivelson et 
al , [1996b]; i.e., the measured field disturbance is wider and 
deeper The model with reflective boundary conditions yields 
a somewhat deeper central disturbance than the fixed boundary 
model, but the overall widths are nearly identical. 

Figure lb shows the plasma density compared with the 
measurements of Frank et al. [1996], Both the peak ot the 
density at the center of the wake and the width ot the 
disturbance are well reproduced by both mass-loading models. 
In fact, the height of the plasma density was one key factor in 
adjusting the mass-loading rate in our process of reproducing 
the measurements. The plasma pressure from the model is 
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compared with the Frank et al. data in Figure lc and the 
temperature in Figure Id. The character of the pressure/ 
temperature disturbance is reasonably well reproduced. There 
is a rise to either side of a narrow central minimum. The depth 
and width of the central minima in pressure and temperature, 
which correspond well with narrow density peak, are \yell 
reproduced by both models. However, the measurements show 
a large increased temperature and pressure In the flanks of the 
wake, which roughly corresponds to the wider measured 
magnetic field disturbance, and which is not quantitatively 
reproduced by either mass-loading model. The models show 
smaller peaks which are closer to the central line of the wake, 
and the model with reflecting boundary conditions produces 
higher side peaks in both temperature and pressure. 

The two key attributes on which we focused fitting were the 
height of the density peak (as mentioned above) and the central 
depth of the minima in temperature and pressure. Although 
there is some interplay in adjusting both parameters, the 
density peak was more sensitive to changes in the mass- 
loading rate, while the depths of the pressure and temperature 
minima were more sensitive to the charge exchange rate. For 


the model results shown in Plates 1 and 2 and Figure 1, the 
total ion r lass-loading rate corresponds to 0.84 x 10 28 and 1.8 
x 10 28 ion> s’ 1 for the fixed and reflective boundary condition 
models, respectively, where the mean molecular mass is 22 amu 
and the radial distribution varied as r 3 5 . The charge exchange 
rate (due again to momentum/energy loading for which charge 
exchange Is the dominant part) is equivalent to 1.5 x 10 28 and 
3.2 x 10 28 ions s’ 1 , respectively, for comparison with the 
charge exchange rates given by Linker et al. [1991]. For the 
best match to the density peak and pressure reversal, the ion- 
neutral collision frequency was about two times the ionization 
rate. The 2 to 1 ratio in relative rates are in rough agreement 
with recent estimates by Smyth [1992]. 

This net ion mass loading corresponds to a fresh ion mass- 
loading rate of 300-650 kg s’ 1 . From the momentum/energy 
exchange rate we can generate an equivalent total apparent 
mass-loading rate which has contributions from both the true 
ionization rate as well as the momentum and energy loading 
from the combination of (mostly) symmetric charge exchange 
and non-reactive ion-neutral collisions of about 840-1800 kg 
s’ 1 . This total loading rate is a factor of 2-4 larger than Smyth's 




Figure 1. Comparison of the mass-loading models with Galileo particle and field measurements. The magnetic 
field measurements in (a) from Kiveison et al. [1996] are shown with the model values. The (b) plasma density, 
(c) pressure, and (d) temperature from the model are shown with thos; from Frank et al. [ 1 996] obtained from 
the plasma ion measurements. The solid lines give the results fortho model with mass-loading and reflective 
boundary conditions. The dashed lines give the results for the mod:! with mass-loading and fixed boundary 
conditions. 
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Mass Loading - Reflective Boundary 




Plate 1 . Magnetic field and streamlines for non-mass-loading and mass-loading flows past Io. (a)-(c) The 
plane aligned along the upstream magnetic field direction. The draped field lines (generally vertical lines) and 
velocity stream lines (generally horizontal lines with arrows) are shown, (d)-(f) The plane perpendicular to 
the upstream magnetic field direction. As indicated. Plates la and Id correspond to the case with no mass- 
loading and reflective boundary conditions. Plates lb and le correspond to mass loading and reflective 
boundary conditions, and Plates lc and If correspond to the mass loading and fixed boundary conditions. 
The spectral color table indicates the magnitude of magnetic field. The strong influence of the Alfvdn wings is 
apparent on the non-mass- loaded flow (Plates la and Id), and is weakened with the addition of mass loading. 
The increased magnetic field perturbation in the wake is apparent in the mass-loaded cases. The spatial scale 
covers from -10 to +12 Io radii horizontally and to ±8 Io radii vertically. 
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No Mass Loading - Reflective Boundary 



Mass Loading - Reflective Boundary 



Mass Loading • Fixed Boundary 



Plate 2. Plasma density for the non-mass-loading and mass- 
loading flows past Io. (a) The plasma density for the models 
with no mass-loading and reflective boundary conditions, in 

(b) for mass-loading and reflective boundary conditions, and in 

(c) for mass-loading and fixed boundary conditions. The non- 
mass-loading model (Plate 2a) produces a rarefied region in the 
wake, contrary to the Galileo particle measurements. Both 
mass-loading models (Plate 2b and 2c) produce a strong 
density peak in the wake behind Io. The symmetry line 
through Io covers distances from -4 to +8 Io radii. 


Voyager-era estimates but may be consistent with the apparent 
overall larger Io atmosphere and torus activity in recent times 
compared vith the Voyager era. 

The general flow pattern in both mass-loading models is 
similar to the typically suggested qualitative pictures. For 
example, see Schneider et al. [1991], Very close to Io and 
especially in the extended wake region the flow is very slow (a 
few km s' 1 ) where the effects from the locally produced ions and 
from charge exchange are most important. Because the Galileo 
trajectory is nearly in the perpendicular plane and because 
there is not much difference between the two models in this 
plane a comparison of the velocity field with the Galileo 
results provides little discrimination between the two 
alternative boundary conditions. Figure 2 shows the flow 
velocity vectors along the spacecraft trajectory for both mass- 
loading models, which give virtually identical results. The 
flow directions and change in the velocity magnitude in the 
orbit plane are quite similar to those shown in Figure 1 of 
Frank et al. [1996], Note again, however, that Frank et al. 
reported asymptotic velocity magnitudes far from the wake (4 
Io radii) somewhat below strict corotation (45 km s' 1 as 
opposed to 56 km s' 1 ). They warned that their uncertainties 
were asymmetric and large enough to include strict corotation, 
being dependent on the details of their assumed plasma ion 
composition. Our model results yield corotational velocities 
far from Io with somewhat supercorotational velocities (13% 
enhancem< nts) in the immediate flanks of the wake near Io, 
about 2 Io radii from the center of the wake. 

In ideal MHD the current density can be calculated as the 
curl of the magnetic field vector. Extracting quantities 
determined from components of field gradients near Io, where 
they are largest, is uncertain because of the presence of cut- 
cells (partial cubes which contain the inner boundary at one 
face), beca lse of the conversion from cel I- averaged parameters 
to node-va ues, and because of boundary effects. For the model 
runs showi : here the cell structure was not optimized to obtain 
highly accurate currents. We can say at least qualitatively that 
the largest values of the current density in all three models 
occur near the body. As one might expect, currents associated 
with the Alfvdn wings are the strongest in the non-mass- 
loading reflective boundary model and weakest for the mass- 
loading fixed boundary model. In the mass-loading fixed 
boundary nodel we can also estimate the integrated current 
from the ci rrent density near Io to be on the order of a few x 10 6 
amperes. This is comparable to independently estimated values 
[Khurana et al. , 1997] for both the Galileo and Voyager 
epochs. 

Preliminary pictures of the global shape of the local 
ionosphere of Io have been presented based on radio 
occupation measurements by Flasar et al. [1997] and plasma 
wave measurements by Gurnett et al. [1997], This overall 
shape is dominated by the dense, stagnant, and heavily mass- 
loaded region indicated by the very low velocities which 
excludes th e outer plasma torus flow, as seen in Plate 2. It is 
apparent hen that the overall structure of Io's upper 
ionosphere may be largely controlled by the combination of 
mass loadi ig and plasma torus flow past Io. Furthermore, the 
upstream-d iwn and down stream -dusk ionosphere profile 

asymmetry seen in the Pioneer 10 radio occultation 
measureme its [Kliore et al , 1974, 1975], that are qualitatively 
similar to he preliminary Galileo measurement, is then most 
likely an upstream/downstream phenomenon rather than a 
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atmosphere and plasma. Including all the effects of this kind of 
detailed picture for the Io flow problems has already been done 
for the comet solar-wind flow problem in terms of successful 
comparisons with in situ spacecraft data [Gombosi et ai , 
1996] and 2-D ground-based spectra and images for H 2 0 + 
[Hdberli et ai, 1997], In the comet problem, although much 
more complicated in terms of total species and reactions than 
for Io, the details of the neutral atmosphere structure and 
composition are also much better understood than for Io. 

The second possibility would be to reexamine the premise 
behind the formulation of the pick-up process. As mentioned 
in the model description, the ion pick-up process adopted 
assumes ions are produced initially in a pick-up ring and 
shortly thereafter scattered into a shell distribution by wave- 
particle interactions. Recent Hubble Space Telescope images 
of the Jupiter aurora with the new STIS instrument [Clarke et 
ai, 1997] show that the footprint of Io is accompanied by a 
trailing footprint of the wake itself extending up to one-quarter 
to one-half way around Jupiter. This implies that ions picked 
up at Io, are still being pitch angle scattered into the loss cone 
for 3 to 6 hours, and still precipitating along wake-crossing 
flux tubes into Jupiter’s ionosphere for a quarter to a half a 
Jupiter rotation. Those picked up very close to Io, where 
densities are high, would be scattered immediately. This could 
also explain why the models provide an excellent match to all 
the measurements in the densest regions. The inclusion of 
pickup ions remaining in ring distributions would 
substantially alter the formulation of the MHD equations. A 
multispecies model that could account for a mixed ring and 
shell populations is beyond our current capability. It is 
certainly well beyond any modeling that has yet been 
published. It is possible that the inclusion of this process (W. 
Smyth, private communication, 1997) might have desired effect 
on the magnetic field disturbance. 


4. Summary 

We have demonstrated the application of an ideal MAUS- 
MHD model to the interaction of Jupiter’s corotating plasma 
torus with Io’s atmosphere. The model solves the single-fluid 
MHD flow in 3-D on an adaptively refined mesh. Because of 
its numerical robustness and the ability to adapt the grid 
structure to the flow, the heavily mass- momentum-energy 
loaded flow in the supersonic and sub-Alfv6nic regime 
relevant for the interaction of Jupiter’s plasma torus has been 
solved without the need to introduce artificial dissipative 
terms. 

We find that we can reproduce many of the observed 
structures of the plasma wake behind Io as measured by various 
Galileo instruments. We obtain a good quantitative match to 
the plasma density peak, the plasma pressure and temperature 
minima at the center of the wake, and the vector velocity plasma 
flow directions .and magnitude variations. We obtain a 
qualitative match to the overall pressure and temperature 
profiles. We reproduce the magnetic field disturbance about as 
well as the intrinsic dipole magnetic field model of Kivelson et 
ai [1996b]. The distant flanks of the disturbance indicate a 
more decreased magnetic field as well as corresponding 
increased plasma pressure and temperature which appear to be 
coincident and possible correlated. Although the responsible 
mechanism is not yet known, it may be due to decreased 


cooling h the flanks from velocity-dependent collision cross 
sections (which are not yet taken into account) or from the 
influence of pickup ions that are still in ring distributions. 
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